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A THIN- SHOCK- LAYER SOLUTION FOR NONEQUILIBRIUM, 

INVISCID HYPERSONIC FLOWS IN EARTH, MARTIAN, 

AND VENUSIAN ATMOSPHERES 

By William L. Grose 
Langley Research Center 

SUMMARY 

* 

An approximate inverse solution is presented for the nonequilibrium flow in the 
inviscid shock layer about a vehicle in hypersonic flight. The method is based upon a 
thin- shock- layer approximation and has the advantage of being applicable to both sub- 
sonic and supersonic regions of the shock layer. The relative simplicity of the method 
makes it ideally suited for programing on a digital computer with a significant reduction 
in storage capacity and computing time required by other more exact methods. 

Comparison of nonequilibrium solutions for an air mixture obtained by the present 
method is made with solutions obtained by two other methods. Additional cases are pre- 
sented for entry of spherical nose cones into representative Venusian and Martian 
atmospheres. 

A digital computer program written in FORTRAN language is presented that per- 
mits an arbitrary gas mixture to be employed in the solution. The effects of vibration, 
dissociation, recombination, electronic excitation, and ionization are included in the 
program. 


INTRODUCTION 

The feasibility of planetary exploration has been the impetus for extensive research 
into problems associated with atmospheric entry of a blunted-nose vehicle. For the 
extreme velocities encountered, the gas flow in the shock layer about the vehicle can 
experience a severe departure from thermodynamic equilibrium. Analysis of the invis- 
cid regime of the shock layer becomes rather complex and requires the simultaneous 
solution of the basic fluid-dynamic equations plus a system of coupled nonlinear equations 
governing the nonequilibrium rate processes. A further complication arises since the 
character of the system of fluid-dynamic partial differential equations changes from ellip- 
tic to hyperbolic as the flow changes from subsonic to supersonic. 


Consequently, numerous techniques, each having its own peculiar advantages and 
disadvantages, have been used to investigate nonequilibrium blunt-body flows. The fol- 
lowing techniques are discussed in reference 1: Newtonian approximations (refs. 2 and 
3); streamtube method (refs. 4 to 6); inverse marching integration (refs. 7 to 10); inte- 
gral relations (refs. 11 to 14); artificial viscosity and time asymptotic finite differences 
(refs. 15 and 16); series truncation (ref. 17); and boundary -layer -type approximations 
(refs. 18 to 20). 

The present method is based upon the thin- shock- layer assumptions of Maslen 
* (ref. 21). It requires specification of a shock shape and solves for the resultant body. 

This disadvantage is inherent in inverse methods, since, in general, the flow over a spec- 
ified body is desired and the shock shape is not known. However, subsequent experience 
demonstrated that many desired body shapes could be obtained with some minor experi- 
mentation with the assumed shock shape. The analysis depends upon casting the governing 
flow equations in shock-oriented coordinates and then performing a Von Mises transfor- 
mation so that the independent variables become the stream function and distance along 
the shock. An approximation is made to the momentum equation in the direction normal 
to the shock that enables it to be readily integrated. From this equation the pressure 
distribution in the Von Mises plane can be determined independently of the chemistry 
within the shock layer. The chemical and vibrational rate equations can then be inte- 
grated along streamlines. Transformation from the Von Mises plane to physical space 
is accomplished by numerical quadrature. 

This approximate method of solution possesses three distinct advantages: 

(1) Applicability throughout the subsonic and supersonic regions of an inviscid 
shock layer; 

(2) Capability of handling a realistic model of an arbitrary gas mixture with a large 
number of chemical species and reactions in terms of reasonable computer memory stor- 
age; and 

(3) Suitability for programing on a digital computer with significant reductions in 
computing times as compared with most existing methods. 

Results obtained by use of this method for nonequilibrium air are compared with results 
from two other computational methods. Additional results are presented for entry of 
spherical nose cones into representative Venusian and Martian atmospheres. 

This solution has been programed in FORTRAN IV and is presented with instructions 
for usage in appendix A by Barbara L. Weigel of the Langley Research Center. The pro- 
gram has the capability of handling up to 50 reactions and 25 species and includes the 
effects of vibrational relaxation, dissociation, vibration-dissociation coupling, electronic 
excitation, and ionization. 
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SYMBOLS 


A coefficient in equation (56) 

Aj frequency factor in rate constant for jth reaction 

a speed of sound, cm/sec 

ay correlation integers for ith species in jth reaction (see appendix B) 

Bj temperature exponent in rate constant for jth reaction 

Cj mass fraction of ith species 

Ej activation energy in rate constant for jth reaction, K 

e internal energy per unit mass, ergs/g 

e g j electronic energy per unit mass of ith species, ergs/g 

ej internal energy per unit mass of ith species, ergs/g 

e v,i vibrational energy per unit mass of ith species, ergs/g 

e v i equilibrium vibrational energy per unit mass of ith species, ergs/g 

fj correlation factor: fj = 0 for monatomic species; f^ = 1 for all other 

species 

g. } degeneracy of Ith electronic energy level of ith species 

AHj standard heat of formation at 0 K of ith species, ergs/mole 

h static enthalpy per unit mass, ergs/g 

h = h/vj 

Kj rate constant for jth reaction, (moles)^" n (cm)^ n "^sec"^, where n is order 

of reaction 
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number of electronic energy levels for ith species 
M Mach number; general species in appendix B 

Mf Mf = 0, vibrational equilibrium; Mf = 1, vibrational nonequilibrium 

nij code indicating which species to use when calculating coupling factor for 

jth reaction 

Nj number of vibrational energy levels included in dissociation energy of 

,, ith species 

n order of reaction 

p pressure, dynes/cm2 


P - P/PeJoo 2 

q^ general parameter representing either vibrational energy per unit mass or 

mass fraction of kth species 

R universal gas constant, 8.3147 x 10? ergs/mole-K 

R c shock radius of curvature, cm 

Rj rate of reaction j per unit mass of mixture 

Rj ratio of forward rate to reverse rate of reaction j per unit mass of mixture 

R n nose radius of body, cm 

r, perpendicular distance from symmetry axis, cm 


r = 


= r/ R N 


moles per unit volume of ith species, moles/cm3 

a* jCj 

S T = 1 if , . = 0; S T = — ^ — if v T i . = 1 

J 1+1,3 ’ J Pi I+lj 
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T 


temperature, K 


T = T/T 00 

T y ^ vibrational temperature of ith species, K 

t time, sec 

u velocity component in x-direction, cm/sec 

V velocity, cm/sec 

v velocity component in y-direction, cm/sec 

x distance along shock, cm 

x = x/R n 

y distance normal to shock, cm 

y = y/y b 

z distance along symmetry axis, cm 

z = z/R n 

ot i jj coefficient in vibrational relaxation time expression 

jj temperature exponent in vibrational relaxation time expression 

dissociation energy of ith species, K 
6 function defined in equation (49) 

e i } energy of fth electronic level of ith species, K 

characteristic vibrational temperature, K 
9 local shock inclination angle (see fig. 1) 
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A 

X 

P 

Mi 



Vl,j 

P 


function defined in equation (47) 

scale factor (1 - y)/R c 

molecular weight of gas mixture, g/mole 

molecular weight of ith species, g/mole 

stoichiometric coefficient of ith reactant in jth reaction 

stoichiometric coefficient of ith product in jth reaction 

stoichiometric coefficient of general species (see appendix B) in jth reaction 
mass density, g/cm3 


P - P/P TO 

a i k coefficient in vibrational relaxation time expression 

T i k vibrational relaxation time of ith species due to vibrational energy exchange 

with kth species, sec 

( p angular position coordinate (fig. 1) 

(p j coupling factor defined in equation (40) 

\p stream function, g/sec 

fl function defined in equation (48) 

*» 

t 

(x)^ • general parameter in rate equation representing either the rate of production 

1 or rate of change of vibrational energy of species k 


Subscripts: 


b body 

s shock 
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oo 


free stream 


0 behind shock on symmetry axis 


ANALYSIS 


Governing Equations 

The solution described herein is an inverse method designed to take advantage of 
the characteristically thin shock layer that forms about a smooth axisymmetric body in 
hypersonic flight. If the shock- oriented coordinate system shown in figure 1 is adopted, 
the steady flow of an adiabatic, inviscid gas subject to nonequilibrium rate processes is 
governed by the following equations: 

Continuity equation: 


l_(pur) + ^(pvAr) = 0 


a 


Momentum equations: 


u 


8u . ... 9u uv 1 8 P 

ax av R„ o ax 


( 1 ) 

( 2 ) 


U |Y +Xv £V + u^ + 2.IK = 0 


av , u 


a 9 P 


3x 


9 y R„ p 9y 


Energy equation: 


h + 



v 2 

— = Constant 


(3) 


(4) 


Rate equation: 

% 9q k 

u ^r +xv iT =Xa, k( p ’ p ’ q i’ • • *’%) 


Equation of state: 


(5) 


h — h^p,p,q^, • • ., (6) 

In view of the thin shock layer, it is convenient to perform a Von Mises transfor- 
mation. The stream function can be defined from the continuity equation by its partial 
derivatives, 

aip 

-r- = pvrx (7) 
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and 


d\p 

w 


-pur 


which leads to the transformation relationship 


( 8 ) 


(d_\ _ d_ XV 

9x U 3y 

The resulting expressions for the momentum and rate equations are 


(: 


8u\ v , Xvr /8p\ _l_/£p\ _ q 

9x/^/ R c u \8^/ x PU\9x/^ 


(9) 


( 10 ) 


\3x/^ R c pj 


(ID 


Xa, k 

\3x Jty u 


( 12 ) 


A differential form of the energy equation results from differentiating equation (4) and 
combining it with equations (2) and (3) 



(13) 


The concept of a thin shock layer now makes certain simplifying approximations 
feasible. The flow within the inviscid shock layer is nearly parallel to the shock itself, 

except near the stagnation point. The term was neglected and the assumption was 

made that X * 1. To this order of approximation, equations (11) and (12) become 





(14) 


aq^\ = 0^ 
MJy u 


(15) 


Adopting the approximation of Maslen (ref. 21), by using the values of r and u imme- 
diately behind the shock, allows equation (14) to be integrated along a normal to the shock 
and results in 


p(x,i//) = P g (x) + 


u s ( x ) 

R c (x) r s (x) 



(16) 
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At the shock 


ib = 


PJS co r s 


(17) 


and on the body the boundary condition is 

i^ b = 0 (18) 

The effect of the approximations is to uncouple equation (11) from the remaining equations 
so that the resulting expression for the pressure (eq. (16)) is independent of the shock- 
layer chemistry and is a function only of the geometry and chemistry of the shock itself. 


From the requirement of constant energy, 


h + 




(19) 


Since the flow is nearly parallel to the shock, the approximation 

v2 - Vq 2 « u2 (20) 

is made, and the following expression for the velocity u results: 

u = [2(h 0 -h)] 1/2 (21) 

In summary, equations (6), (13), (15), (16), and (21) form a determinate system of equa- 
tions that can be integrated stepwise along streamlines to yield p, p, h, u, and q k 
in the Von Mises plane. In addition to the boundary conditions of equations (17) and (18), 
the oblique shock relations (ref. 22) must be satisfied. 

By noting the geometry of figure 1, the transformation back into physical (x,y) space 
can be accomplished by integrating equation (8) at constant x. The resulting expression 



can be readily evaluated numerically. 


( 22 ) 


The most important implication of the preceding analysis is that the equations are 
now applicable uniformly from the subsonic region through the supersonic region, except 
along the symmetry axis where the transformation is not well behaved. 


Thermodynamic Model for a Reacting Gas 

The gas under consideration is assumed to be a mixture of ideal gases. (See 
ref. 23.) The internal energy of the gas includes contributions from translation, rotation, 
vibration, dissociation, electronic excitation, and ionization. 
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The thermal equation of state is given by 


with the mixture molecular weight given by 



The mass fraction Cj is subject to the constraint 

I 



i=l 

The internal energy of the ith species can be expressed as 


( 23 ) 


(24) 


(25) 



where the rotational (rigid rotor) and electronic energy modes are assumed to be in equi- 
librium with the translational mode. In the computer program to be described subse- 
quently, two options exist for computing the vibrational (linear harmonic oscillator) 
energy. In the first option, the vibrational mode is assumed to be in equilibrium with 
the translational mode and e v j is calculated from 


e 


v,i 


fiRGi 


Mi 



- 1 


(27) 


The second option is for vibrational nonequilibrium and e v ^ is determined by integra- 
tion of the vibrational rate equation along a streamline. 

Since a mixture of ideal gases has been assumed, the total energy is 

I 

e = ^ (28) 

i=l 
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The enthalpy is defined as 


h = e + £ (29) 

By referring to the previous equations, it can be seen that 

h = h(p;p;e v i , . . e v>I ;c lf . . Cj) (30) 

which agrees with the form stated in equation (6). 


Rate Equations 

The general rate expression giv^n in equation (5) must be specialized for chemical 
processes and for vibrational relaxation. The resultant equation must be capable of 
describing a system of chemical species (i = 1,2,3, . . ., I) and the associated kinetic 
reactions (j = 1,2,3, . . ., J). 

Chemical processes .- An arbitrary one-step chemical reaction can be described by 

I 

£ (31) 

i=l 


X 

I 

i=l 


v i S i 


The law of mass action (ref. 23) requires the rate of production of a species to be propor- 
tional to the product of the concentrations of the reacting species, each concentration 
being raised to a power equal to its stoichiometric coefficient. The net rate of produc- 
tion of per unit mass of mixture is then 




(32) 


Since 



(33) 


equation (32) can be generalized for a system of J simultaneous reactions and put in the 
form of equation (12). As a result, 



(34) 
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The total rate of production is 



( 35 ) 


The specific reaction rate constant is given by 


Bi 


/-E 


Kj=AjT Jexp^- 


(36) 


Equation (34) is applicable to opposing reactions if the forward and reverse reactions are 
treated separately. 


If the rate constants for the reaction of a chemical species with several catalytic 
species are the same, provision can be made to incorporate all these catalysts into a 
single reaction by introducing a "general species." (See appendix B.) The resulting 
rate expression is now 



(37) 


Vibrational relaxation .- If the vibrational rate equation developed in reference 24 
for a multicomponent gas mixture is written in the form of equation (12), the result is 



(38) 


where e y ^ is the equilibrium vibrational energy from equation (27) and e v> i is the 
vibrational energy. The relaxation time T. . can be expressed as 

1,K 


T i,k~ 


Pi k 
a i,k T ’ ex P 


(<u*- ,/3 ) 


i r e i 
1 - exp — i 


(39) 


which is of the form predicted by classical Landau-Teller theory. Although equations (38) 
and (39) were formulated on the assumption of a system of harmonic oscillators having a 
single vibrational mode, polyatomic molecules with several vibrational modes can be han- 
dled in the computer program if it is assumed that they are independent of one another 
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and the relaxation frequencies for the individual modes are combined into a single equiv- 
alent frequency. (See ref. 23.) 

Coupling of vibrational relaxation and dissociation .- It has been pointed out in the 
literature (for example, refs. 25 to 28) that the processes of vibrational relaxation and 
dissociation do not occur independently, but rather each process influences the other. 

The CVD (coupled-vibration-dissociation) model of Hammerling et al. (ref. 25) was 
chosen for this analysis. Without discussing the relative merits of the various other 
coupling models, it should be noted that another model could be incorporated into this 
method with little difficulty. 

Since molecules with vibrational energy less than that corresponding to equilibrium 
with the translational energy do not dissociate as readily as predicted by the dissociation 
rate constant for vibrational equilibrium, the CVD model develops a correction to the rate 
constant. From reference 25 the coupling coefficient is 


1 - exp 


Vjjp 






expta^l - 1 
v,i. 


( e i\ i 

exp (■=*) - 1 


(40) 


The number of vibrational states Nj is the smallest integer so that 
A.- 

Ni > -i 

x ©i 


(41) 


Inserting the coupling coefficient in equation (37) results in 



For a dissociation reaction, (p j can be calculated from equation (40); for other reactions, 
( p . is unity. 


SOLUTION PROCEDURE 
Initial Conditions 

An inverse method presupposes specification of the shock shape. In all cases 
described herein, an analytical expression of the form 
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= r s (z) 


( 43 ) 


was used. 


From this relationship, 


Q, R c , and x can be given by 


9 = arc tan 



and 


R c - 


r 

1 + 



3/2 


d2 r s 

dz2 


(44) 


(45) 


=i; 


S dr s 
sin 9 


(46) 


After determining the shock parameters, initiation of the solution begins with cal- 
culation of free- stream conditions. With selected values of p T V and a sneci 

fied gas composition, equations (23) to (28) can be utilized to calculate a p and e 

The next step requires determining the state of the gas immediately behind the 
shock wave. The shock wave is assumed to be a discontinuity, translational, rotational 
and electronic contributions to the internal energy of the gas reaching instantaneous equi- 
1 num at the translational temperature behind the shock T s . The vibrational energy 
may be treated as being in equilibrium with the other energy modes or optionally it may 
be described by a system of harmonic oscillators with a Boltzmann distribution about a 

nonequilibrium vibrational temperature. Species concentrations are assumed to be frozen 
through the shock. 

An iterative procedure is required to solve the oblique-shock-conservation rela- 
tions From the shock geometry and the free-stream conditions, the following parameters 
can be calculated : 


A(x) = p V sin 9 

CO 

&(*) = Poo + P oo V oo 2 sin 2 0 




'00 p 


(47) 

(48) 

(49) 


With an initial estimate of T s , the energy e s can be determined from equations (26) to 
(28). The temperature and energy are further related by 
14 


( 50 ) 


2 P<x >( 6 ■ e s)| 

/ft 2 - 2 A 2 (6 - 

e s) 

R 

ft + ^ft 2 

1 - 2A 2 (6 - e s ) 



Equation (50) can be iterated until subsequent values of T s differ by a sufficiently small 
amount. After convergence 


P s = - 2 A 2 (6 - e s ) 

= P S f i o 0 
P s RT S 


h s 


+ 


P s 


(51) 

(52) 

(53) 


u s = cos e 


(54) 


The pressure distribution in the x,i// plane can now be calculated from 
equation (16). 


Integration Method 

Since the initial conditions and the pressure distribution have been determined, 
equations (13), (37), and (38) can be integrated stepwise along streamlines. It was found 
that for the cases investigated herein, the system of equations exhibited the so-called 
"stiff" behavior. - (See refs. 29 and 30.) The integration scheme proposed by Treanor 
(ref. 31) was chosen for its ability to handle the stiff equations with significant reductions 
in computational times compared with both the fourth-order Runge-Kutta method (ref. 32) 
and the Adams-Moulton method using the fourth-order Runge-Kutta method for starting 
(ref. 33). 

A special problem occurs for the \f/ = 0 or body streamline. On all other stream- 
lines the integration begins immediately behind the shock and proceeds from there. How- 
ever, on the Z-axis the transformation of the governing equations is ill-behaved since x 
and \p are both zero. Therefore, the integration cannot proceed from behind the shock 
along the stagnation streamline with the equations in their present form. 

The i// = 0 streamline integration was accomplished in the following manner. The 
initial interval Ax on the shock wave was divided into five subintervals. The rate equa- 
tions were integrated from behind the shock at each subinterval to the point x = Ax. The 
values of the variables as a function of if/ were then extrapolated along x = Ax to the 
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body (i// = 0). Integration along the \]s = 0 streamline can proceed from there. Five 
subintervals were decided to be adequate after experimenting with up to 15 subintervals. 

Transformation to Physical Space 

Upon completion of the integration along streamlines, the thermodynamic properties, 
the velocity u, and the species concentrations are known functions of x and \p. Equa- 
tion (22) can then be evaluated numerically by Simpson's rule to obtain r. Evaluation of 
the flow field is thus completed. 

RESULTS AND DISCUSSION 

In order to investigate the capabilities of the present method of solution, cases have 
been computed for numerous shock shapes and gas mixtures. A number of representative 
cases have been selected for presentation. Case I is a comparison of results for nonequi- 
librium air with the more exact inverse method of reference 8. Case n is a comparison 
of the results for nonequilibrium air with the streamtube method used in reference 34. 
Cases in, IV, and V present results for entry of spherical nose cones into representative 
Martian and Venusian atmospheres. 


Case I 

This example was computed to compare with the results of reference 8 for a cate- 
nary shock (R c> 0 = 2-11 cm) with a free-stream velocity of 7.01 x 10^ cm/sec at 61 km. 
The reactions and rates given in appendix C are the same as those used in reference 8. 
Vibrational equilibrium was assumed. 

Figure 2 indicates excellent agreement with the results of reference 8 for the 
body shape and the location of a streamline crossing the shock at an angle 6 of 81° 
(x/R c q = 0.16). The variations of the temperature, pressure, density, and species mole 
fractions along this streamline are compared with the corresponding results from ref- 
erence 8 in figures 3, 4, 5, and 6, respectively. As should be expected, the effect of 
neglecting the velocity gradient term in equation (11) is to predict the pressure to be 
lower than the actual value at a given x. This effect can be expected to diminish as 
the flow proceeds away from the stagnation region and becomes more nearly parallel 
to the shock. 

Initially, the flow along this streamline exhibits behavior similar to that behind a 
normal- shock wave with density increasing with sharply falling temperature. However, 
as the flow expands away from the nose, the density reaches a maximum and begins to 
decrease as the temperature continues to decrease. 
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Concentrations of N, O, NO, and e~ shown in figure 6 indicate very good agreement 
with the results of reference 8. The greatest error occurs for NO shown in figure 6(c), 
the present method predicting the concentrations about 20 percent too high at worst. 

Variations of the pressure, temperature, density, and species concentrations across 
the shock layer at 9 = 81° are shown in figures 7, 8, 9, and 10, respectively. Again, the 
agreement with the results of reference 8 is good. From this comparison it may be con- 
cluded that the present method gives good agreement with the exact method near the stag- 
nation region which represents a severe test of the approximations inherent in the present 
method. 


Case n 

This example was selected for comparison with the nonequilibrium air results of 
reference 34 for flow over a 9° half-angle cone with a spherical nose (30.48-cm diameter) 
at an altitude of 47.5 km. The equation 

r g 2 = -0.3625z g 2 + 36.064z s - 0.1767r s z s (55) 

was fitted in the nose region to the shock points computed in reference 34. Although Evans, 
and others, used a sophisticated system of 13 species and 54 reactions, the approximate 
air model used in case I was used for this comparison. Vibrational nonequilibrium was 
assumed. 

It can be seen in figure 11 that the present method predicts the body shape extremely 
well. A comparison of the surface pressures and temperatures as calculated by both 
methods is presented in figures 12 and 13, respectively. In general, the agreement is 
seen to be good. There is a discrepancy in the temperature calculations close to the 
symmetry axis, but this condition improves rapidly with increasing distance. Both 
methods are admittedly inaccurate in the region where the discrepancy occurs. Finally, 
the electron concentration on the body surface is shown in figure 14. The prediction of 
the present method is understandably lower since the reaction scheme permits only the 
NO + ion. Agreement improves further back on the surface. This result might be expected 
since figure 7 of reference 34 indicates NO + is the major ion present after the flow expands 
past the sphere- cone tangency point. 


Case HI 

This case was computed for ballistic entry into a Venusian atmosphere, 
found that shock shapes of the form given by 

, -Isi- 

s 1+Ar s 


It was 


( 56 ) 
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would generate bodies that approximate part of a sphere-cone. The parameter A repre- 
sents the slope of the asymptote of the shock for large values of r s . For sufficiently 
large apex angles under appropriate conditions, the flow over the face of a truncated 
sphere-cone is subsonic, the sonic point being on the body at the sharp corner where the 
transition to supersonic flow occurs. There is a resultant rapid drop in pressure near 
the corner. It should be noted that a shock shape of the form of equation (56) can be used 
to calculate the flow in the shock layer up to the point where the rapid pressure drop 
occurs, but another equation would be required for the shock where it "bends" as the 
flow expands around the sonic corner. However, in this rapidly expanding flow region, 
the approximations inherent in this method of analysis would not be valid. 

According to recently available information on Venus (ref. 35), the atmosphere is 
almost entirely carbon dioxide. Therefore, a 100-percent C0 2 mixture was assumed 
and the kinetic reactions and rates shown in appendix C were used. Vibrational nonequi- 
librium was assumed. The free-stream conditions are 

Voo = 8.69 X 105 cm/sec 

P to = 1.4175 x 10-6 g/ cm 3 

T— =-200-K 

The assumed shock (A = 1.732) and the computed body which approximates a 57° half- 
angle sphere- cone are shown in figure 15. Pressure, density, and temperature variations 
along streamlines crossing the shock at x = 0.5 and 2.0 are presented in figures 16, 

17, and 18, respectively. Figure 16 indicates that the pressure along both streamlines 
approaches the same constant level as the flow progresses away from the shock as would 
be expected for cone flow. 

Species concentrations along these streamlines are shown in figure 19. Figure 19(a) 
indicates that the C0 2 disappears rapidly with increasing distance while the C and O con- 
centrations increase along the inner streamline (x = 0.5). The CO concentration increases 
to a maximum and then decreases. Figure 19(b) shows similar behavior along the outer 
streamline (x = 2.0). This behavior is the same as that observed by Howe, Viegas, and 
Sheaffer (ref. 36) behind normal shocks in C0 2 at slightly different ambient conditions! 
Their gas model however did not include the effects of vibrational relaxation or ioniza- 
tion. Figure 19(c) indicates that the concentrations of the electrons and the C+ and 0+ ions 
increase with distance whereas CO + reaches a maximum and then begins to decrease along 
the inner streamline (x = 0.5). The same behavior is exhibited along the outer streamline 
(x = 2.0) except that the CO + has not reached a maximum and begun to decrease, as shown 
in figure 19(d). Along both streamlines, the dominant electron- producing reaction is the 
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ionization of CO. The C 2 concentrations shown in figures 19(e) and 19(f) are several 
decades lower than the concentrations of the other species. 

The pressure and temperature distributions along the surface of the body are shown 
in figures 20 and 21, respectively. Figure 22 presents species concentrations along the 
body surface. On the body surface the longer particle residence times have allowed more 
CO dissociation and resulted in most of the electrons being formed from ionization of C 
and O. This result is in contrast to the results shown along the inner and outer 
streamlines. 

Figure 23 is a profile of the static enthalpy across the shock layer at x = 2.0. It 
can be seen that the enthalpy increases sharply near the body. This region near the body 
is known as the entropy layer. Gas particles in this region have experienced a much 
larger rise in entropy from traversing the shock closer to the stagnation streamline and, 
as a result, much of the kinetic energy was converted to thermal energy. Pressure and 
temperature profiles across the shock layer at x = 2.0 are shown in figure 24. A typi- 
cal electron concentration profile across the shock layer is shown in figure 25. It can be 
seen that the concentration builds up rapidly within the entropy layer near the body. 


Case IV 

This example was selected to represent entry into a possible Martian atmosphere 
having a composition of 95 percent CO 2 and 5 percent A by mass. The free-stream con- 
ditions are: 

V M = 4.572 x IQS cm/sec 


p = 1.068 X 10-5 g/cm3 

OO 

Too = 200 K 


The scheme of kinetic reactions and rates used in this example is given in appendix C. 
Vibrational nonequilibrium was assumed. After some experimentation it was found that 
a shock shape of the form 


r 2 
x s 

Zs “ 1 + 1.89r s 


(57) 


resulted in the body shown in figure 26 which closely approximates a 60° half-angle 
sphere-cone with a 1-cm nose radius. Pressure and temperature variations along 
streamlines crossing the shock at x = 0.4 and 1.6 are shown in figures 27 and 28. 


Species concentrations along these same streamlines are presented in figure 29. 
The CO 2 concentration decreases rapidly with increasing distance as it dissociates into 
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CO and O and results in the increasing concentrations of these species shown in fig- 
ure 29(a). It is seen from figure 29(c) that the dominant ion present is CO + . 

Figure 30 illustrates the pressure and temperature variations along the surface of 
the body. Species concentrations along the body are presented in figure 31. It can be 
seen that the peak electron concentration on the surface occurs near the nose. The elec- 
trons are seen to result primarily from ionization of CO as was the case on the stream- 
lines within the shock layer. 

The sudden rise in static enthalpy within the entropy layer is shown in figure 32 for 
x = 1.6. Pressure and temperature profiles across the shock layer at this value of x 
are shown in figure 33. Figure 34 illustrates the variation in the electron concentration 
across the shock layer. A marked increase in electron concentration is observed within 
the entropy layer. 

Investigation of the assumed reaction mechanism for this case indicates that the 
chemical kinetics are controlled predominantly by two-body collision processes. Since 
the three -body recombination reactions are relatively unimportant, the remainder of the 
reactions are two-body collisions having the same density dependence. For specified ini- 
tial temperature and species concentrations, the time (or distance) variation of all vari- 
ables can then be scaled with respect to density. Practically, this concept of binary 
scaling (ref. 8) means that for a given flight velocity, flows will be similar if the product 
of ambient density and nose radius is held constant. For this particular case the results 
would then apply for a body with nose radius of 30.48 cm and p^ = 3.5 x 10“ ^ g/cm3. 
These conditions might be typical of a Viking- type entry on Mars. 

To demonstrate the validity of the binary scaling assumption, it is necessary to 
assess the importance of three-body recombination reactions for this system. Define the 
rate of reaction j per unit mass of mixture as 



and define the ratio 

R j,j+1 = ^/^j+l ( 59 ) 

The ratio Rjj+i can then be used to judge the relative importance of the forward rate 
to the reverse rate of a reaction. 

Shown in figure 35 are the ratios along the inner streamline (x = 0.4) of the forward 
rate for the system of reactions as numbered in appendix C for this case. It can be seen 
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that Rj j + i » 1 for all the reactions involving three-body recombinations except R3 4 
which indicates that the recombination of C and O to form CO is becoming more impor- 
tant far from the shock on this streamline. However, at the furthermost point, R3 4 ~ 6; 
therefore, the assumption that the flow along this streamline can be scaled is probably 
still a reasonable approximation. A true test of the validity of binary scaling in this case 
would require repeating the calculation for various values of p M for constant p K R n 
and comparing the results. 


Case V 


This example is similar to Case IV with the same free -stream conditions and kinetic 
reactions, but with a composition of 50 percent CO2 and 50 percent A by mass. After 
some experimentation it was found that the shock shape given by 


r s 2 


Zs 1 + 1.98r s 


(60) 


results in the body shown in figure 36 which closely approximates a 60° half-angle sphere- 
cone with a 1-cm nose radius. The shock standoff distance is greater for this case than 
that for Case IV because of the lower density ratio across the shock for this mixture. 
Pressure and temperature variations along streamlines crossing the shock at x = 0.4 
and 1.6 are shown in figures 37 and 38. Comparing figures 28 and 38 shows that the tem- 
perature in the shock layer is about 10 percent higher for this case. The variations of 
species concentrations along these streamlines are presented in figure 39. Comparison 
with figure 29 from Case IV shows essentially the same behavior. Again, the predomi- 
nant ion is CO + . The electron concentrations are somewhat higher for this case because 
of the higher temperature in the shock layer. 

Pressure and temperature variations on the body surface are shown in figure 40. 
Figure 41 illustrates variations of species concentrations on the surface. The only obvi- 
ous difference in the profiles from those of Case IV appears in the electron concentration. 
In the previous case the profile reached a maximum near the sphere- cone juncture and 
then decreased. For this case the profile decreases monotonically with increasing 
distance. 


Existence of the entropy layer is again demonstrated by the rapid increase in static 
enthalpy near the surface shown in figure 42 for x = 1.6. Pressure and temperature dis- 
tributions across the shock layer at this value of x are shown in figure 43. The varia- 
tion, in electron concentration across the shock layer at this value of x is shown in fig- 
ure 44 to be about a factor of 10 higher than that for Case IV. 

The ratio of forward reaction rates to reverse reaction rates for this reaction sys- 
tem as given in appendix C is shown in figure 45. It can be seen that all Rj » 1 for 
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the three body recombination reactions. Thus, the binary- scaling assumption appears to 
be justifiable for this case. 


CONCLUDING REMARKS 

An approximate inverse solution is presented for the nonequilibrium flow in the 
inviscid shock layer about a vehicle in hypersonic flight. A number of representative 
cases have been presented. Case I is a comparison of results for nonequilibrium air 
with the more exact inverse method. Case n is a comparison of the results for nonequi- 
librium air with the streamtube method. Cases m, IV, and V present results for entry 
of spherical nose cones into representative Martian and Venusian atmospheres. From 
the foregoing analysis and results, the following observations can be made. 

The method of solution as programed herein can handle a realistic model of an arbi- 
trary gas mixture with large numbers of chemical reactions and species. A larger num- 
ber of reactions and/or species could be easily accommodated by increasing the dimen- 
sioning of the appropriate variables in the FORTRAN program. The gas model permits 
consideration of vibrational relaxation, dissociation, recombination, ionization, electronic 
excitation, and vibration-dissociation coupling. With relative ease the gas model could 
be made more sophisticated. For example, the vibrational partition function could include 
effects of anharmonicity or a more realistic vibration-dissociation coupling model could 
be included. 

Many flow-field solutions suffer the disadvantage of being valid only in the subsonic 
(or supersonic) region of a shock layer. Another method is required to solve the super- 
sonic (or subsonic) region so that the solutions match in the transonic zone. The approxi- 
mations inherent in the method presented herein avoid this difficulty and render the 
method applicable for subsonic and supersonic flows. 

The solution offers little difficulty in programing for use on a digital computer as 
evidenced by the program described herein. Primarily because of the approximations 
made to the governing equations, this method offers the possibility of significant reduc- 
tions in computing time. The required computing time for Case I on the Control Data 6600 
digital computer was 79 seconds. The remaining 4 cases required between 5 and 14 min- 
utes of time, Case n requiring the maximum time. These times are based upon calcu- 
lating 20 streamlines. This result appears to be a significant reduction in the times 
reported in most of the existing literature on nonequilibrium flow fields. 

The accuracy of the present method compared extremely well with the inverse 
marching integration technique for Case I and the streamtube method for Case n. 

Although the method fails on the Z-axis and the approximations are crude close to the 
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stagnation point, the requirements of a thin shock layer and a smooth symmetric body 
are not so restrictive that application of the method is unduly hindered. 

For Case in the electron concentrations were found to be greatest near the body in 
the entropy layer. In this region the electrons are produced primarily from ionization of 
C and O atoms. Within the shock layer, but outside the entropy layer, the concentrations 
are lower and the major ion present is CO + . 

For Case IV the electron concentrations were found to be greatest in the entropy 
layer. Throughout the shock layer the major ion present is CO + . Maximum electron 
concentration occurred near the nose on the body. 

The assumption of binary scaling for this system of reactions used in this case 
appeared to be a justifiable approximation. It was evident, however, that the recombina- 
tion of C and O to form CO is the critical reaction in determining the validity of this 
assumption. 

The larger percentage of argon for Case V compared with Case IV results in the 
argon acting as a nondissociating diluent with resulting higher temperatures in the shock 
layer. The electron concentrations are correspondingly higher. Maximum concentra- 
tions exist in the entropy layer. The binary -scaling assumption again appears to be jus- 
tifiable for case V. 

It should be emphasized that the systems of reactions chosen for these cases are 
merely representative and by no means definitive. Possible inaccuracies may arise if 
an incomplete reaction mechanism is used. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va., October 28, 1971. 
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APPENDIX A 


A COMPUTER PROGRAM TO DETERMINE THE FLOW FIELD IN 
A NONEQUILIBRIUM SHOCK LAYER ABOUT A SMOOTH 
SYMMETRIC BODY TRAVELING AT HYPERSONIC 
SPEEDS USING AN INVERSE TECHNIQUE 

By Barbara L. Weigel 
Langley Research Center 

Introduction 

This appendix describes the digital computer program (LRC Program D1290) devel- 
oped to support the study of an inverse technique (shock- shape specified) for determining 
the flow field in a nonequilibrium shock layer about a smooth symmetric body traveling at 
hypersonic speeds. The problem description, development of governing equations, sym- 
bols, method of solution, and results are discussed in the text of this paper. 

Problem Description 

The problem of determining the flow field in a shock layer about a smooth symmet- 
ric body behind a specified shock shape may be divided into six parts : 

(1) Computation of the shock geometry (x s , z s , r s , R c , cos 6, and sin 0) from a 

definition supplied by the user. (See fig. 1) 

(2) Computation of the free-stream quantities P> a <»> and 

(3) Computation of quantities immediately behind the shock (e s , T s , p g , p g , u g , 

and \p s ^ 

(4) Computation of pressures, and their derivatives, in the flow field (p (x,i//) and 

(5) Integration along streamlines from a point just behind the shock to a point suffi- 

ciently downstream to define the flow field 

(6) Computation of physical space quantities r, y, and z, at specified points in 

the flow field to give a good body and flow field definition 
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APPENDIX A - Continued 
Numerical Methods 

The numerical methods used for the development of this program include 

(1) A modified Runge-Kutta integration scheme designed for chemical problems 

(ref. 31) 

(2) The Newton- Raphson iteration method (ref. 37, p. 192) 

(3) Numerical differential equations (ref. 38, p. 96) 

(4) Simpson's Rule integration technique (ref. 39, p. 137) 

(5) Trapezoidal Rule integration technique (ref. 39, p. 142) 

Subprograms 

The subprograms used are given in the following listing: 


FORTRAN 

name 

Called 

by 

Function 

SHOCKG 

D1290 

to compute the shock geometry and 
store on tape 10 

FDPDX 

D1290 

to evaluate the pressure p and its 
derivative p and store on tapes 

FOFTS 

ITR1 

to evaluate Tg 

ITR1 

BASIC 

to iterate for e 

ITR1 

D1290 

to iterate for T s 

CALINTH 

D1290 

to obtain the integration technique 

BASIC 

CALINTH 

to evaluate the derivatives h, c i( and 
e v i for i = 1,1 

FOFE 

ITR1 

to evaluate e 

CHECK 

CALINTH 

to test results of an integration step and 
accept or reject the results 
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APPENDIX A - Continued 

Subprogram SHOCKG(DELX,ZSTERM,IZTERM) is called by D1290 to compute the 
shock geometry parameters and store them on tape 10. This subroutine must be supplied 
by the user to define the shock shape desired. The calling sequence must contain the two 
input quantities DELX and ZSTERM and a generated integer, IZTERM, the number of x 
values computed in the subroutine. The user must store the following parameters on 
tape 10: x, z, r, R c , cos 9 , and sin 0. These six parameters must be computed at 
each x from x = 0 to x ^ ZSTERM and x must be in increments of DELX/5 to 
DELX and increments of DELX thereafter until x = ZSTERM where 

x distance along shock, cm 

z distance along symmetry axis, cm 

r perpendicular distance from symmetry axis to shock, cm 

R c radius of curvature of shock, cm 

9 local shock inclination angle, radians 

STOP 13 has been used for error stops in this subroutine. In the sample case, the fol- 
lowing shock shape was defined: 

z = cosh r - 1.0 


x = \/2z + z2 
cos 9 = x 

\ PTlo 

R c = 1.0 + x2 

Iteration of e 

To encourage iteration to an acceptable value of e and thereby reduce instability, 
the following tests have been used: 

(1) If the maximum number of iterations (10) is exceeded before convergence is 
reached, a new estimate of e is computed and the iteration is again attempted. The 
program will restart this iteration twice. If convergence is not reached, the program 
will terminate with STOP 66. 

(2) If the converged e is greater than h and a negative T results, e is set 
equal to h - 0.000001 * h and one more iteration on e is attempted. 
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APPENDIX A - Continued 


(3) If, after this maneuver, e is again greater than h, a trigger is set to restart 
the integration interval at one-fourth the current computing interval. 


Testing Results of Integration 

After an integration interval has been computed, several tests are made to decide 
whether to 

(1) Continue, 

(2) Recompute at one-fourth the current computing interval, or 

(3) Stop the case. 

Trial and error has proved that recomputing at a smaller interval sometimes prevents 
instabilities and error conditions from developing. 


Stagnation Streamline 

The stagnation xfy = 0.0 streamline is computed after the streamlines 



2 Ax 3 Ax orlH 4 Ax 
— — and — 


have been integrated to x = Ax. Values of e, T, h, Cj, and e ± (i = 1,1) at xf/ = 0.0 
are obtained by extrapolation. Then T v i (i = 1,1), u, p, p, and p values are com- 
puted and the i p = 0.0 streamline is integrated from Ax to x = ZSTERM. 

Negative u2 

In the evaluation of the derivatives if u2 becomes negative, the integration is dis- 
continued and an attempt is made to compute the physical space values for the calculated 
flow field. 


Input 

The data input is loaded by using FORTRAN IV NAMELIST. One card, 80 columns, 
is read in by FORMAT (8A10) for case identification. The input symbols are as follows: 


FORTRAN 

name 

$NAM1 ' 

DELX 

ZSTERM 


Symbol 


Description 


Ax increment along shock, cm 

length of symmetry axis, z, cm 


27 


APPENDIX A - Continued 


FORTRAN 

name 

Symbol 

Description 

IMAX 

I 

maximum number of species i, §25 

JMAX 

J 

maximum number of reactions j, §50 

MJ 

m i 

an array containing an integer indicating 
which species i to use to calculate the 
coupling factor <p j for reaction j 

M 

Mf 

1, vibrational nonequilibrium; 0, vibra- 
tional equilibrium 

R 

R 

universal gas constant, ergs/mole- K 

GAMMA 

y 

ratio of frozen specific heats 

CIINF 

c i>°° 

an array containing free- stream mass 
fraction for each species i 

PINF 

Poo 

free- stream pressure, dynes/cm2 

TINF 

Tco 

free- stream temperature, K 

VINF 

Vqo 

free- stream velocity, cm/sec 

MUI 

fH 

an array containing molecular weight for 
each species i, g/mole 

BJ 

B i 

an array containing temperature exponent 
for each reaction j in Arrhenius -type 
rate equation 

EJ 

E i 

an array of activation energies for each 
reaction j , K 

AU 

a i,j 

a two-dimensional array of correlation 
factors for ith species in jth reaction 
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FORTRAN 


name 

Symbol 

Description 

NUU 


a two-dimensional array of stoichiometric 
coefficients of ith reactant in jth reaction 

NUPU 

ij 

a two-dimensional array of stoichiometric 
coefficients for ith product in jth reaction 

ALPIK 

“i,k 

a two-dimensional array of coefficients in 
vibrational relaxation time expression 

BETAIK 

%,k 

a two-dimensional array of temperature 
exponents in vibrational relaxation time 
expression 

SIGK 

a i,k 

a two-dimensional array of coefficients in 
vibrational relaxation time expression 

THETAI 

®i 

an array containing characteristic vibra- 
tional temperature K for each species i 

DGENI 

Bi 

an array of factors for each species i to 
permit approximating a polyatomic mole- 
cule by a diatomic molecule 

FI 

*i 

an array of correlation factors for each 
species i, 0.0 for monatomic species 
and 1.0 for all others 

DELHI 

A Hi 

an array of standard heat of formation at 
0 K for each species i, ergs/mole 

DELI 

A i 

an array of dissociation energy for each 
species i, K 

EVI 

e v,i 

* an array of vibrational energy per unit 
mass for each species i, ergs/g 
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FORTRAN 

name 

LI 

GIL 

EPSIIL 

AJ 

XI 

ELE1 

ELE2 

XPST 

NXPST 

EPF 


Symbol Description 

Lj an array containing number of electronic 

levels for each species i 

gj j a two-dimensional array containing degen- 

eracy of LI electronic levels for each 
species i 

e. , a two-dimensional array containing energy 

of LI electronic levels for each species 
i, °K 

Aj an array containing frequency factor, 

cm^ n /mole n -sec, for each reaction j 
in Arrhenius -type rate equation 

the initial computing interval. If no XI 
is input, 0.0001220703125 will be used 

an array of (2(IMAX) + l] values used by 
integration scheme to control size of 
computing interval. 0.0 < ELE1 ^ 64.0 

an array of (2(IMAX) + l| values used by 
integration scheme to control size of 
computing interval. ELE2 < ELE1 

an array of x values at which physical 
space calculations are desired. They 
must be multiples of DELX. The program 
sets XPST(l) = DELX and XPST(NXPST) 

= gx at ZSTERM) + 100. 6] 

number of values in XPST array = 100 

PRINT ANSWERS every IPF integration; 

IPF = 5 unless input otherwise 
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FORTRAN 

name 

Symbol 

Description 

CIMAX 


maximum computing interval; if CIMAX is 
not input, 0.0625 will be used 

TIMER 


time, in decimal seconds, requested on 
JOB card. If case does not terminate 
in time estimated on JOB card, data for 
pickup will be stored on tape 13. If 
TIMER is not input, 3600.0 or 1 hour 
will be used 

IPICKUP 


= 1 to pickup case from tape 13 
= 0 unless input 

IBUG 


an array of 15 integers which must be 


input as 0 to trigger debug printouts 
at critical points (see Program Listing 


$ 

With some cases, difficulties were encountered in starting the integration. When 
this condition occurred, the input quantities XI, ELE1, ELE2, and CIMAX were varied 
until the integration proceeded. 

After reading in NAMELIST input, read in the IMAX species (10 spaces per species, 
8 species per card) by Format (8A10) and the JMAX reactions by the same FORMAT using 
20 spaces per reaction, 4 reactions per card. If the species and reactions are not read 
in, IMAX/8 + JMAX/4 blank cards must be put in the deck. If (IMAX/8) and (JMAX/4) 
are not whole numbers, they should be rounded to the next whole number. 


Sample Input 


A sample input follows: 

SHOCK SHAPE NEQA GROSE ♦ 12- 17-68 7 SPECIES. 14 REACTIONS OOOO 

SNAM1 
w*bt 

MU I (I )= 28* 0 1 6 1 32 • OOO ♦ 14*000* 16*000« 30 • 008 * 30* Q08 • 5 •♦862E— 4 • 

Cl INF (1 >=0#768« .232, .0 • • 0 ♦ *0 • .0 * • 0 ♦ 

PINF=1 #98E2* 

TINF * 253*8« 

VINF = 7*02E5* 
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R =8 
GAW«A= 1 
THETA I ( 1 
GIL(1 *1 
GIL< 1 *2 
G I L ( 1 .3 
GIL ( 1 *4 
G I L ( 1 ,5 
GILC l 1 6 
G I L ( 1 1 7 


DGENI (1 >® 20*1 *0 


EPS! IL( 1 
EPSI IL( 1 
EPSI IL( 1 
EPSlILd 
EPSI IL( 1 
EPS ML C 1 
EPSI IL C 1 
LI <1 > = 4 
FI ( 1 ) = l 


31 469E7. 

4. 

a 3392# 14, 2275# 0*1 *0.1 .0.2739.0.3421 #0,1 #0* 


* 1 *0 « 3 

* 3.0,2 

38 4 • .6. 

* 5.. 3. 

= 2 .. 2 . 

3 1 . . 6 . 


= 1 • ♦ 7 * 0 • 


1 )* 0 . 
2 ) = 0 . 
3 ) = 0 * 
4 ) * 0. 

5 ) * 0 . 

6 ) ® 0 . 

7)® 8*' 

6 * 4 • 5 * 6 

. 1 • .0# . 


AJ( M 
BJC1 > = 
B J (11 
EJ(1 > = 
E J (11 
NUI J C 1 
NUI J( 1 
NUI Jt 1 
NUI JC 1 
NUI J( 1 
nuijT 1 
NUI J( 1 
NUI J( 1 
NUI J( 1 
NUI J( 1 
NUI JC ! 


A I J ( I 
A I J ( 1 
A I J ( 1 
A I J ( 1 
A I J ( I 
A I J ( 1 
AI JC 1 
A I J C 1 


-“1 *5* — 1 .0.-1 #5.- 
* - 1 • 5 • 0 . 1 2 .0 ♦ 0 . 

59380.. O.. 7549 

® 0 • ♦ 3 1858 • *0 * 0 . 


0 « 6 • 0 • 2 • 0 ♦ 4*0 • • 

0.1 .0 .3.0. 1 .0. 3. 0.2*0 • • 

4 . . 6. .4*0* ♦ 

1.. 5..1. .3*0 . ♦ 

2. . 4. .2. .4. .2. .4.. 

3 • .5*0 . . 


71592 .. 85343 .. 99212 # . 4 * 0 ## 

1 1342 .. 18879 # . 51 385 # . 52104 # # 71026 # • 2 * 0 *. 
27659 #. 27670 # . 41496 #. 4 * 0 *. 

228 .. 326 #. 22831 #. 48622 *. 3 * 0 #. 

174 *. 63597 • . 6538 1 • . 76676 • • 7536 1 * • 87568 • * 86359 • • 
57528 #. 842 05 #. 5 * 0 # ♦ 


3.1. 

#••1 # . 1 • *0. . 

DELHI Cl )«6#.0#.4#707E12.2.4681EI2.9#0073E1 1 . 1 • 0826E 1 3 . 0 • . 

DELI Cl)® 1 14981 #*60500#*0.#0#. 76883# .127698# .0*. 

M J ( 1 ) ® 2. 4. 5.3. 1.3*2. 5. 1.5*1 * 5 .6.3# 

MJC 1 1 >« 6* 3.0.0. 

AJC1 )= 1 #2E21 . 1 .0E18.5.2E21 .1 #3E21 «3#0E21 .1 #67E20. 1 #0E1 2.2 #38E 1 1 «5#0E1 3. 
1 .1 1E13.9.1E24.4.79E23. 1 #8E21 .2.17E1 1 • 


1 .1 1E13.9.1E24.4.79E2: 
I® 1 #8E21 .2.17E11 .0.0. 


— 1 #5.-1 #5.-1 #5. — 1 #5.0#5.0#5.0# .0# .**2 #5 2#5# — l #5*0# 12. 


1 )= 0.1.0 
2 ) = 0.0.0 
3 ) ® 0.0.0 

4 ) ® 0*0*1 

5) ® 1.0.0 

6 ) = 0. 6T2 

7>= 0.1.1 
8)® 0.0.0 
9 ) ® 1*0.0 

to )= 0 . 0.1 

11)= 0.0.0 


75490. #0*. 1 1 3260. .0# . 3120#. 19130#. 3600O#.0#« 

. 0 . 0 . 


0 • 0. 0 , 0 * 1 . 
2.0. 0,0*1 « 
0.1 .0,0*1 . 

1 . 0 . 0 . 0*1. 
0.0. 0,0 * 1 . 

0 . 0 . 0 . 0 * 1 « 
0 , 0 , 0 , 0 * 0 . 
1 «r« o,o«o. 

1 .0. 0. 0*0. 
0.1 .0,0*0. 
0.0.1«1«0. 


NUI JC 

1 * 

1 

2)® 


0. 

0. 

1 . 

1 ,0,0, 0*0, 

NUP1 J 

C 1 


1 ) 

a 

0 

.0 

.0 

.2.0. 0.0. 

NUP! J 

Cl 


2) 

a 

0 

.1 

♦ 0 

. 0.0.0, 0. 

NUPI J 

c 1 


3) 

a 

0 

.0 

. 1 

• I .0.0.0. 

NUP! J 

c 1 


4) 

a 

0 

.0 

.0 

♦0.1 .0.0. 

NUP! J 

(1 


5) 

a 

0 

,0 

.2 

•0 .0.0.0. 

NUPI J 

c 1 


6) 

a 

1 

♦ 0 

.0 

.0.0.0 .0. 

NUP I J 

c 1 


7) 

a 

0 

♦ 0 

.0 

.1.1 .0.0. 

NUP I J 

(1 


8) 

a 

0 

,1 

. 1 

• 0 .0 .0.0. 

NUPI J 

c 1 


9) 

a 

0 

.0 

. 1 

•0.1 .0*0. 

NUPI J 

< l 


10) 

a 

1 

.0 

.0 

. I .0.0.0 • 

NUPI J 

C 1 


11 ) 

a 

0 

.0 

• 1 

.1.0. 0*0. 

NUPI J 

c 1 


12 ) 

= 

0 

• 0 

.0 

.0.0.1 ,1 . 

A I JC 1 

• 

1 

> = 

6*1 

.0 

. 


AI JC1 

• 

2) = 

6*1 

.0 

. 



3 ) = 6*1 

4) t* 6*1 

5) = 6*1 


1 0 ) = 7*0 
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A! J(1 .11 )* 7*0, 

A!J(1» 12 ) = 7*0 , 

ALPtKtli t)= 7*5.675E-4. 
ALPIKd. 2>* 7*1.443E-7. 
ALPIKd. 3 ) * 7*0, 

ALPIKd « 4)* 7*0, 
ALPIKth 3)* 7* I • 51 3E— 4 « 
ALP IK (!• 6)° 7*1 .51 3E— 4 . 
ALPI K (1, 7>° 7*0. f 
BETA IK ( 1 ♦ 1 ) *7*-. 09952 f 
BETA I K < 1 • 2 ) “7*0 • 371 77 f 
BETA IK ( 1 f 3 ) *7*0 . • 

BETA I K ( 1 . 4 ) *7*0 • * 

BETA I K ( 1 . 5>*7*-.2226, 
BETA ! K d . 6)»7» -.2226 . 

BETA IK U , 7) *7*0, 
SIGIKClf l )= 7*163.63* 

S IGIK < 1 , 2> a 7*179.32. 
SIGIKC1 , 3 ) 3 7*0. 

S IGI K ( 1 • 4> a 7*0. 
SIGlKdf 5)= 7*135.3. 

S IGI K ( 1 • 6>= 7*135.3. 
SIGlKdf 7 ) * 7*0, 

IMAX=7. 

JWAX= 1 2 • 


CIMAX=.5. 

ELEI d 1*42*64., 

ELE2 C 1 ) *42*1 .0, 

X I * .0625, 

ZSTE9W»1 .145. 
0ELX=0.06. TCHECKT*) 
EVI < 1 ) =20*0. 0 , 

NXPST *15. XPST C 1 ) *= • 06 
I PF*5 . fIWER*240.f 
S 

M2 02 

02 + M 3 20 + M 

N2 + M ■ 2N + M 

N2 + O = NO + N 


.3. HCHECKT* 0 #01 

.12 , * 18. .24. .3. .36 
IPICKUP=0, 

N O 

20 + M = 02 + M 

2N + M a N2 + M 

MO + N * N2 + 0 


, PHMAX=65*0* 

. 42 • .48 » . 54 , , 6 • • 66 . 


MO NO+ 

no + m*n + o-+m 

N + 02 = NO + O 
N0+ + E- = N + 0 


.72. .78. .84.10.. 


E- 

N + O+ MaNO + M 
MO + O a N + 02 
N + O * NO* ♦ E- 


Tape Assignments 
Tape 5, input, with buffer size 1001 
Tape 6, output, with buffer size 1001 

Tape 8, to save derivative of pressure at each x for each \p 
Tape 9, to save pressure at each x for each i// 

Tape 10, to save the shock geometry parameters x and corresponding z, r, 
R c , cos 0, and sin 0 

Tape 11, to save quantities behind the shock T s , e s , p g , p s , u g , \f/ s , h g , 

( e v i) s ^ = at each x 

Tape 12, to save physical space data 


Tape 13, pickup tape 
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Operational Details 

Program D1290 was written in FORTRAN IV language for the Control Data 6000 
series digital computer under the Scope 3.0 operating system. The program requires a 
field length of 110 000g. The sample case was completed in two passes on the Control 
Data 6600 computer to illustrate the pickup capability. 

Programed STOPS 

STOP 1, incorrect input 

STOP 6, EOF 10 in main program 

STOP 13, in SHOCKG 

STOP 30, in CHECK when computing interval less than 1.0E - 15 
STOP 50, EOF 9 in FDPDX 

STOP 66, in BASIC when no convergence on e iteration 

STOP 301, in MAIN for error in XPST array or IZTERM .GT. 500 

STOP 321, ITR1 stop in main program 

STOP 422, in MAIN after pickup tape written 

STOP 663, in MAIN when x is not equal to VARI(IPSI) 

STOP 665, in MAIN after an integration attempt 
STOP 670, when ac^ is negative 

Pickup Tape 

Under control of the input value TIMER, the estimated time on the JOB card and a 
REQUEST card, a pickup tape will be generated if the case does not terminate in the esti 
mated time. A pickup may not be made before the = 0.0 streamline integration has 
started. 

Output 

When no debug printouts are requested (IBUG array = 1), the printed output will 
contain: 

1 ID information 

2 a list of the IMAX species 

3 a list of the JMAX reactions 
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4 a printout of the input quantities and those generated in the initialization 

section 

5 the shock geometry x, z, r, R c , cos 6, and sin d 

6 the free-stream quantities | u^, p^, a w) e ioo , e M , ( e e,i) » and 

(e v>i j^ (i = 1,JMAX) 

7 the quantities behind the shock T s , e s , p g , p g , u s , i// g , h s , and 

(e v (i = 1,IMAX) for each streamline or each x distance 
along shock 

8 IZTERM, the number of x values on the shock or the number of streamlines 

to be integrated 

9 The following values are printed every IPF acceptable integration intervals : 


IPSI 

streamline ID beginning with 2 for xj/ at x = DELX/5 

IX 

number of integration steps for streamline beginning with 1 at 


the shock 

x,h,p ,p^ 
P»u,T ,ej 

IMAX 

r — ' 

for computing interval / c^ 

i=l 

c t (i=l,EMAX) 

concentration 


c m,i (i=l,IMAX) molecular concentration 
e ,, : (i=l,IMAX) vibrational energy 

V 9 1 

(IPSI=6) is used for the stagnation streamline, x(/ - 0.0, and the streamline begin- 
ning at Ax 

10 after all the streamlines have been computed, the physical space values r,y,z 

are printed at the x values given in the XPST array 

11 after the initial conditions for a streamline are printed, the statement 

" CHECKT= " appears. This number represents the elapsed 


time, in seconds, from start 
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12 other comments are printed to help the user follow the flow of the solution or 

to explain built-in remedial maneuvers to continue a solution which is in 
trouble 

13 the IBUG debug printout options represent areas of difficulty during the 

development of the program. They may be helpful to the user for a closer 
examination of selected quantities and so have been left in the program 

Sample Output 

A sample output follows: 
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APPENDIX A - Continued 


X- 6.00000000E-01 
x= 6.00000000E-01 
X® 6 • 0000G000E-01 
X® 6 *0000 0000 E -01 
X- 6*00000000 E-01 
X= 6 • 0000 0000 E-01 
X® 6 .OOCOOCOOE-OI 


X= 6.60000000E-01 
X= 6.60000000E-01 
X= 6 # 6000 0000 E -01 
X= 6 • 6000CC00E -01 
x= 6.60000000E-01 
X= 6 *60000000 E -01 
X- 6.60000000E-01 
X- 6 • 60000000E-0 1 
X® 6 * 600 00 000 E -0 1 
X® 6*60000 OOOE-0 1 
X- 6.6000C 000E-01 
X- 6.60000000E-01 


X® 7 • 20000000E-01 
X- 7*20000000E-01 
X® 7*200000 00 E-01 
X- 7*200000006-01 
X® 7.20000000E-01 
X® 7.20000000E-01 
X= 7* 20000 OOOE-0 1 
X= 7.2OCOC0OOE-O1 
X= 7 • 2000C000E-01 
x= 7 • 20000000E-01 
X= 7 • 2000 C OOOE-0 1 
X® 7.20000000E-01 
X® 7*20000000E-01 


x= 7*900000 00 E-01 
x= 7* 80000000 E-01 
X= 7 • 80 000 OOOE-O 1 
X® 7*8000CC00E-0i 

x= 7 *80000000 E-01 

X® 7 • 8 000 0000 E-01 
X® 7* 800 000C0 E-01 
x= 7.80000000E-01 
X- 7.80000000E-01 
X= 7 • 80000 OOOE-O 1 
X® 7 • 8 0000 OOOE-O l 
X- 7*8000 00 00 E-01 
X= 7 • 80000000E-01 
X» 7 * 8 0000 OOOE-O 1 


X- 9*400000 00 E-01 
x= 8 *40000 000 E-01 
X= 8*4CC OCOCOE-Ol 
X- 8.40000000E-01 
X® 8 • 400 00 OOOE-O 1 
X* 9*400000 00 E-01 
X= 8 *400000008 -Cl 
_X = 8*40000 000E-01 
X=” 8 • 40000000 E-0 1 
X® 8*400000008-01 
X® 8*400000008-01 
X= 8 .40000000E— 01 
X® 8. 400000 00 E-01 
X® 8*400000 00 E-01 
X= 8* 400000 OOE-Ol 


END THIS CASE, 


R = 5* 57248692 E-01 
R - 5* 61 75571 7E-0 1 
R® 5*668498046-01 
R® 5*724608016—01 
R* 5.78566507E-01 
R = 5*851182 40 E-01 
R® 5*919702996-01 


R = 5. 9909653 3 E-01 
P® 5* 996 662 8 5 E-01 
R® 6*01289080 E-01 
R * 6*038843436-01 
R® 6* 0736764 4 E-01 
R * 6*11608163 E— 01 
R® 6. 16422309E-01 
R« 6*21761643 E-01 
R® 6*276039026-01 
R® 6* 3391 871 9E-01 
R® 6* 4062 790 3 E-01 
R = 6*471475076-01 


R® 6.490088726-01 
R * 6*49536208E-01 
R® 6*51050279 E-01 
R= 6*534854746-01 
R® 6*567657906-01 
R® 6.607758116-01 
R® 6.6534611 IE-01 
R® 6*704357656-01 
R® 6*7602817 5E- 01 
R= 6. 8 20 97234 E-01 
R® 6*885968046-01 
R® 6*949391356-01 
R® 7*013040836-01 


R® 6*98778666 E— 01 
R= 6*992715946—01 
R= 7*006985826-01 
R® 7 *0300441 6 E-0 1 
R® 7*0612 34 88 E-01 
R= 7.099479556-01 
R® 7* 14322554E-01 
R® 7.19210973E-01 
R® 7. 24591602 E-01 
R® 7* 3044 9774 E-01 
ft® 7*367519326-01 
R® 7* 42949722 E-01 
R® 7.49120716E— 01 
R= 7.5546C289E-01 


R® 7*483086296-01 
R® 7.487745936-01 
R® 7*501329176-01 
R= 7*523346506—01 
R® 7 • 55321601E-01 
R® 7*58997210 E-01 
R= 7*63217648 E-01 
R= 7* 67945558 E-OJ. 
R® r.T3i 59275E-0 1 
R= 7*788439936-01 
R® 7*8497 64 11 E— 01 
R= 7 • 91 041936E-0 1 
P® 7*9712 744 8E-0 1 
P= 8*032740856—01 
R= 8* 09616659E-01 


Y= 1 . 2702 203 7E-01 
Y® 1 • 10531504E-01 
Y= 9*1893 002 8E-02 
Y= 7 • 1363201 5E-02 
Y® 4.90233357E-02 
Y= 2*50’51 5239E-02 
Y®-1 *913549876-05 


Y® 1 *625439596-01 
Y= 1*606164446-01 
Y® 1.551264436-01 
Y= 1 • 4634651 6E- 01 
Y® 1*345 62306E-01 
Y= 1 • 20216380E-01 
Y= 1 *0392 983 9E -01 
Y= 8 • 5 866 556 3E -02 
Y® 6*6101 85 13E-02 
Y= 4*473 84543 E-02 
Y® 2*204088886-02 
Y=-1.53163277E-05 


Y® 1 *651 6083 8E-01 
Y® 1 . 6349526 1 E-0 1 
Y= 1 • 587 131 02 E-01 
Y® 1 * 51021600E-01 
Y® 1.40660800E-01 
Y= 1 • 27995247E— 01 
Y= 1 .135600676-01 
Y= 9*748451626-02 
Y= 7 *982 1024 4E -02 
Y® 6*0652 05 13E-02 
Y= 4*0123314 8E- 02 
Y = 2.00912163E-02 
Y®-1.23167902E-05 


Y® 1*678741016-01 
Y= 1 • 6641 41 1 9E— 01 
Y= 1.62187590E-01 
Y= 1*5535 80 4 8E -01 
Y® 1*46119810E-01 
Y® 1.34792296E-01 
Y= 1*218 3 53 74E-01 
Y= 1 . 07356591E-01 
Y= 9.14199524E-C2 
Y = 7 • 4068902 9 E-02 
Y® 5 • 5402 82 94 E-02 
Y= 3* 7045 87 96 E-02 
Y= l*87682977E-02 
Y=- 8*5936 56 27E-06 


Y® 1 * 709 3 80 6 4 E-01 
Y® 1 * 6963 882 8E -01 
Y= 1 *6585 1442 E-01 
Y= 1 • 5971 239 9E -Cl 
Y® 1 • 5 13 8 395 OE -01 
Y® 1*4113 53 29E-01 
Y= 1*2936 75 75E-C1 
Y= 1.16184846E-01 
Y® 1 • 0164755 7E-C1 
Y= 3.57969824E-C2 
Y® 6* 8698096 3E -02 
Y= 5.17857271E-02 
Y= 3* 481 7 62 63E -02 
Y® 1*7 6790947E-02 
Y=-5.76366770E-06 


Z= 2.0 6759283E-01 
1 = 1 *90904276 E-0 1 
Z= 1* 72984081 E-01 
Z® 1.532 45 475 E-01 
Z= 1* 31766562 E-01 
Z= 1*08718596 E-01 
Z= 8.46141310E-02 


Z= 2*61872375 E-01 
Z= 2*60035972 E-01 
Z= 2.54805480E-01 
Z® 2*46440572 E-01 
Z= 2.35213393E-01 
Z® 2.21545588E-01 
Z= 2 *06028900 E-01 
Z= 1*88 8 19454E-01 
Z= 1*69 989012 E-01 
Z= 1*496 35446 E-01 
Z= 1.28010777E-01 
Z= 1*069 97 152 E-01 


Z= 2*880874276-01 
Z= 2*065 16659E-01 
Z® 2*820067116-01 
Z® 2*747530266-01 
Z= 2 • 649 81 985 E-01 
Z® 2*530373 83 E-01 
Z= 2.39423835E-01 
Z= 2.24263389E-01 
Z= 2*07605342 E-01 
Z® 1*995275 09 E-01 
Z® 1 • 7016731 3E-01 
Z® 1 * 51275485E-01 
Z= 1 • 323 16289E-01 


Z® 3* 14341 92 5E-01 
Z= 3. 12979142 E-01 
Z= 3*09033995 E-01 
Z= 3*026591 30E-01 
Z® 2*940359296-01 
Z® 2 *83462 546E-01 
Z® 2.71368233E-0I 
Z= 2.57853377E-01 
Z= 2.42977722E-01 
Z= 2*26781 821 E-01 
Z® 2*093 584 44 E-Cl 
Z® l *92223611 E-01 
Z® 1* 7 5162 861 E-01 
Z® 1*57636044 E-01 


Z® 3 *40 8 68 009 E-01 
Z= 3.396678UE-01 
Z= 3*36169130 E-01 
Z= 3.30498055E-01 
Z® 3.22804468E-01 
Z® 3*133370796-01 
Z= 3*024663596-01 
Z= 2.90288526E-01 
Z® 2. 768593846-01 
Z= 2.62217C69E-01 
Z® 2*464215986-01 
Z= 2.30798427E-01 
Z= 2*151237726-01 
Z® 1.99291677E-C1 
Z= 1*82954902 E-01 


PSI= 5*44413801E-03 
PSI® 8.47881 147E-03 
PSI® 1.21346707E-02 
PSI® 1*64440873E-02 
PSI® 2* 14070611 E-02 
PSI® 2*70235921E-02 
PSI® 3.32936805E-02 


PSI® 0. 

PSI® 3.41887158E-04 
PSI® 1*366591556-03 
PSI® 3.06679722E-03 
PSI® 5* 44413 801 E-03 
PSI® 8*4788U47E-03 
PSI® 1.21346707E-02 
PSI® 1*64440873 E-02 
PSI® 2 • 1407061 IE-02 
PSI® 2* 70235 92 IE-02 
PSI® 3*32936805 E-02 
PSI® 3*978965146-02 


PSI® 0* 

PSI® 3.41887158E-04 
PSI® 1*36659155 E-03 
PSI® 3*066797226-03 
PSI® 5 *4441 3 801 E-03 
'PSI® 8*47881 147 E-03 
PSI® 1*21 34 6707 E-02 
PSI® 1.64440873E-02 
PSI= 2.14070611E-02 
PSI® 2*70235 921E-02 
PSI® 3.32936805E-02 
PSI® 3* 97896514E-02 
PSI® 4*672803866-02 


PSI® 0. 

PSI® 3.418871586-04 
PSI® 1 * 366591 55E-03 
PSI® 3*066797226-03 
PSI® 5* 444138 01 E-03 
PSI® 8*478811476-03 
PSI® 1*213467076-02 
PSI® 1 *64440 873 E-02 
PSI® 2*140706116-02 
PSI® 2.702359216-02 
PSI® 3*32936805 E-02 
PSI® 3.97896514E-02 
PSI® 4.67280386E-02 
PSI® 5*422375536-02 


PSI® 0. 

PSI® 3*41887158E-04 
PSI® 1.36659155E-03 
PSI® 3 *06679722 E-03 
PSI® 5*444138016-03 
PSI® 8.47881147E-03 
PSI® 1.21346707E-02 
PSJ= 1* 644408736-02 
PSI® 2 • 14070611 E-02 
PSI® 2.70235921E-02 
PSI® 3.32936805E-02 
PSI® 3.97896514E-02 
PSI® 4.67280386E-02 
PSI® 5*422375536-02 
PSI® 6.22768016E-02 


HALLELUJAH 
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APPENDIX A - Continued 


Complete Program 

The complete program including comments for the sample case, the flow diagram, 

and the deck setup is reproduced herein. 

PROGRAM D1290 . I INPUT =1001»OUTPUT.=1 001»TAPE5=INPUT ,TAPE6= OUTPUT , 

I T AP E8 » TAPE9, TAP E13 » TAPE11 »T APE1 2, TA PE 13 I 

TAPE8 TO SAVE DERIVATIVE OF PRESSURE AT EACH X FOR EACH PSI 
TAPE9 TO SAVE PRESSURE AT EACH X FOR EACH PSI 
TAPEIO TO SAVE X AND CORRESPONDING ZS. RS.RC, COST, SINT 
TAPE11 TO SAVE T S, ES ,P S ,RHOS ,US, PS I S ,HS . ( E VI S< I > , I = 1 , 1 MA X ) 

TAPE12 TO SAVE PHYSICAL SPACE DATA 
TAPE13 TO SAVE COMMON FOR PICKUP TAPE 

BL WEIGEL FOR WILLIAM L GROSE, APO, GAS PHYSICS SECTION 

INVERSE TECHNIQUE FOR DETERMINING THE FLOW FIELD IN A SHOCK 
LAYER ABOUT A SMOOTH SYMMETRIC BODY TRAVELING AT HYPERSONIC 
SPEEDS 

PHASE 3, A DIRECT NON-EQUILIBRIUM CALCULATION 
CALINTH, INTEGRATION ROUTINE 

THE NUMERICAL INTEGRATION ALGORITHM USED IS FOUND IN -A METHOD 
FOR THE NUMERICAL INTEGRATION OF COUPLED FIRST ORDER 
DIFFERENTIAL EQUATIONS WITH GREATLY DIFFERENT TIME CONSTANTS- 
BY CHARLES E. TREANOR, CONTRACT N0.NASR-119.CAL REPORT NO. 

AG- 1729-A-V. JANUARY 196*. CORNELL AERONAUTICAL L ABORATORY , INC. 

CORNELL UNIVERSITY, BUFFALO, N.Y. 

ITR1 - NEWTON-RAPHSON ITERATION METHOD 
- FTLUP - INTERPOLATION ROUTINE 

R .SHOCK 

• 

1 • 

B SHOCK GEOMETRY 1 • 

IX. • BODY 

1 / . 

1. 

1 • 

• ... . • . . Z 

1 . 

1. 

1 . . 

1 . . 

1 • 

X . 

input" 

READ IN X CARD, 80 COLUMNS, BY FORMAT ( BA10 > FOR CASE ID 
INPUT DEFINITIONS NAMELIST INPUT 


SNAM1 

DELX = INCREMENT ALONG SHOCK 
ZSTERM= LENGTH OF SYMMETRY AXIS ,Z 

IMAX = MAXIMUM NUMBER OF I-S, SPECIES, LESS THAN OR = 25 
JMAX = MAXIMUM NUMBER OF J-S, RE ACT I ONS , LE SS THAN OR = 50 
IBUG = 0 DEBUG PRINT OUTS 

= X , NO DEBUG PRINT OUTS, 1 UNLESS INPUT OTHERWISE 
MJ = CODE INDICATING WHICH SPECIES. I. TO USE_ TO CALCULATE 
COUPLING FACTOR, PHI SUB J, FOR REACTION J 
M = X FOR VIBRATIONAL NON-EQUILIBRIUM. 

= 0 FOR VIBRATIONAL EQUILIBRIUM 
R = UNIVERSAL GAS CONSTANT, ERG/MOLE DEG K 

GAMMA = RATIO OF SPECIFIC HEATS 

Cl INF = FREE STREAM MASS FRACTION FOR EACH SPECIES 
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PINF = FREE STREAM PRESSURE, DYNES/ CM**2 
TINF = FREE STREAM TEMPERATURE, DEG.K 
VINF = FREE STREAM VELOCITY, CM/SEC 
MU I = MOLECULAR WEIGHT FOR EACH SPECIES, GM/MOLE 

THET AI = CHARACTERISTIC VIBRATIONAL TEMPERATURE , DEG. K 
DGENI = FUDGE FACT DR TO PERMIT APPROXIMATING POLYATOMIC MOLECULE 
BY 0 1 ATOM IC MOLECULE 
FI » 0 FOR MONATOMIC SPECIES 
= 1 FOR ALL OTHERS 
DELHI = HEAT OF FORMAT ION, ERGS/MOLE 
DELI = DISSOCIATION ENERGY OF SPECIES 
EVI * VIBRATIONAL ENERGY OF SPECIES 

LI = NUMBER OF ELECTRONIC LEVELS FOR EACH SPECIES LI.LE.20 
GIL = DEGENERACY OF L-TH ELECTRONIC LEVEL FOR I-TH SPECIES 
EPSIIL* L-TH ELECTRONIC ENERGY LEVEL FOR I-TH SPECIES 
AJ = FREQUENCY FACTOR IN ARRHENIUS TYPE RATE EO. 

BJ = TEMPERATURE EXPONENT IN ARRHENIUS TYPE RATE EO. 

EJ = IN ARRHENIUS TYPE RATE EQ. 

A I J * FACTOR TO ALLOW USE OF GENERAL SPECIES IN REACTION EQS IF 

DESIRED, =1.0 OR =0.0 

NUIJ = STOICHIOMETRIC COEFFICIENTS OF I-IH REACTANT IN J-TH 
REACTION 

NUPIJ = STOICHIOMETRIC COEFFICIENT FOR I-TH PRODUCT IN J-TH 
REACTION 

ALPIK = FACTORS IN EQ. FOR VIBRATIONAL RELAXATION TIME 

bTtaTk= fa“ctoFs~Tn - eq. for "vibrational relaxation - ti"me 

SI Gl K .= FACTORS IN EQ. FOR VIBRATIONAL RELAXATION TIME 
XI = INITIAL COMPUTING INTERVAL, .0001220703125 UNLESS INPUT, 
ELE1 = ( 2 * IMAX* 1 1 VALUES USED BY INTEGRATION SCHEME 
NORMALLY .1 ,.5 OR .05 

ELE2 = ( 2* IMAX+1 ) VALUES USED BY INTEGRATION SCHEME 
NORMALLY .05, .1 OR .01 AND.LT.ELE1 
XPST * 99 OR LESS X-S AT WHICH PHYSICAL SPACE CALCULATIONS ARE 
DESIRED. THEY MUST BE MULTIPLES OF OELX IN ORDER TO 
HAVE RS .COST ,ZS AND S INT VALUES 
AND LAST MUST BE .GT.X AT ZSTERM 

THEREFORE, XPSTINXPSTI SET = X AT ZSTERM ♦ 100.0 IN 
PROGRAM 

XPSTIl) MAY NOT BE 0.0 
THEREFORE, XPST(1> SET * OELX IN PROGRAM 
NXPST = NUMBER OF X-S AT WHICH PHYSICAL SPACE CALCULATIONS ARE 
DESIRED 

IPF = PRINT FREQUENCY, 5 UNLESS INPUT OTHERWISE 

IPF NOT IN COMMON THEREFORE IT MAY BE VARIED ON PICKUP RUN 
CIMAX = MAXIMUM CJ OR COMPUTING INTERVAL 
0.0625 UNLESS INPUT OTHERWISE 

CIMAX NOT IN COMMON THEREFORE IT MAY BE VARI EO ON PICKUP 
RUNS 

TIMER = TIME, IN DECIMAL SECONDS, REQUESTED ON JOB CARD 

IF THE JOB DOES NOT TERMINATE IN THE TIME ESTIMATED ON 
THE JOB CARD, DATA FOR PICKUP WILL BE STORED ON TAPE 13. 
IT IS = 3600.0 OR 1 HOUR UNLESS INPUT OTHERWISE 
TIMER NOT IN COMMON THEREFORE IT MAY BE VARIED ON PKKUP 
RUNS 

IPICKUP = 0 UNLESS INPUT 

= 1 TO PICKUP CASE FROM TAPE 13 
HCHECKT = CONTROL ON SIZE OF COMPUTING INTERVAL IN CHECK 
IFIABSIHPR EV-H I / H. GT .HCHEC Kl REDUCE INTERVAL 
TCHECKT = CONTROL ON SIZE OF COMPUTING INTERVAL IN CHECK 
IF ( ABSI TPREV-TI /T.GE. TC HECK ) REDUCE INTERVAL 
PHMAX = CONTROL ON COMPUTING INTERVAL .LE.65.0 
$ 


AFTER READING IN NAM1 , READ IN THE IMAX SPECIES BY FORMATI8A10I 
10 SPACES PER SPECIES 

THEN READ IN THE JMAX REACTIONS BY THE SAME FORMAT USING 20 
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C SPACES PER REACTION 

C IF THE SPECIES AND REACTIONS ARE NOT READ IN, 

C I MAX/B PL uf J MAX/A 8 LA NK CAP 0 S MUST ” BE PuT IN THT~DECiT 

C 

C STOPS 

C STOP 1 INCORRECT INPUT 

C STOP 6 EOF 13 IN MAIN 

C STOP 13 IN SHOCKG 

C STOP 30 IN CHECK WHEN COMPUTING INTERV AL. LT. 1. OE-1 5 

C STOP 50 EOF 9 IN FDPDX 

C STOP 66 IN BASIC WHEN NO CONVERGENCE ON E ITERATION 

C STOP 301 IN MAIN FOR ERROR IN XPST ARRAY OR IZTERM.GT. 500 

C STOP 321 ITRl STOP IN MAIN 

: C STOP A 22 IN MAIN AFTER PICKUP TAPE WRITTEN 

|C STOP 663 IN MAIN WHEN X. NE. V AR I ( IP SI J 

C STOP 665 IN MAIN AFTER AN INTEGRATION ATTEMPT 

C STOP 670 IN MAIN WHEN A Cl NEGATIVE 

C 

COMMON M,IX,IMAX,IPSI,IBUG(15) 

COMMON P< 5001 ,OPOX( 500),VARI( 5001 
C FOLLOWING 7 VARIABLES DIMENSIONED BY I MAX 

COMMON EV 1 1 NF ( 25 ) * THET AI (25 ) ,MUI ( 25 ) ,F I ( 25 ) , DELH I ( 25 I ,CI I NF< 25), 

1 L I ( 25 ) 

COMMON R, MUINF, DELTA, LAMSQ.OMEGSO. OMEGA, MU, T,U 
,C FOLLOWING 2 VARIABLES DIMENSIONED BY (LMAX IN LI,IMAX» 

COMMON EPSIIL(20,25t,GIL(20, 251 
COMMON VAR(52) ,CUVAR(52> ,DER(51> 

C 

COMMON ZS*RS»RC*SINT , COST, PS , US, PS IS 
C FOLLOWING A VARIABLES DIMENSIONED BY JMAX 

COMMON MJ ( 50 ) » EJ (50 ) , AJ ( 50 ) , BJ (50) 

C FOLLOWING 3 VARIABLES DIMENSIONED BY I MAX 

COMMON TVI(25),NI(25>,DGENI(25) 

FOLLOWING 3 VARIABLES DIMENSIONED BY ( IMAX+1 • JMAXI OR 

UMAX * JMAX) 

COMMON NUIJ(26,50),AIJ(25,50) , NUPI J ( 25,50) 

FOLLOWING 3 VARIABLES DIMENSIONED BY (IMAX,IMAX) 

COMMON SIGIK(25,25) , ALP IK( 25 ,25) *BET AIK( 25 , 25 1 
COMMON KSTAG,IMAXP2. IMAXL.N, IG,I STAG,LPS2,LPSI .IMAXP1 
COMMON X,TS,ES,RHOS, HS »RHOURT , PHMAX, AAAAA, PREPSIS,PR6X,PRERU 
EVIS(IMAX) 

COMMON E VIS (25), XPST (100) 

COMMON /LAB2/ DELX ,ZSTERM, IZTERM 

CO MMON / LA B3/ EINF , P INF , RHOI NF, VINF, E, JMAX. KEY INT.RHO. HST AG 
COMMON /LA 84/ PFTL.K ITRl ,E ITRl - 
COMMON /LAB5/ IUNEG 

COMMON /LAB6/ ELB, SPEC »CJ,TPREV,HPREV, HCMECK,TCHECK 
COMMON /LABT/ ITNEG,IEXP 
COMMON /LABS/ EL El (51 ) ,ELE2< 51 ) ♦ NERR 

SPEC IE(IMAX), REACT (JMAX) «OEL I ( IMAX ) 

DIMENSION SPECIE (25 ), REACT (50) , INI 8 ) ,OEL t ( 25) 

DOUBLE PRECISION VAR^XVAR,H,CI ,EVI 

DIMENSION RUT ( 200 ) , PS I STG ( 200 I 
DIMENSION ENG(U),ENGI(11) 

DIMENSION CIGI25.6) , EIG( 25,6) 

01 ME NS ION PSIG ( 5 ) , TG ( 5 ) , HG (6 ) , CIT( 6 ) ,E IT ( 6 ) , EG ( 5 ) 

CM(IMAX),CI( IMAX ),EVI( IMAX) 

DIMENSION CM( 25 ),CI(25 ),EVI(25 ) 

DIMENS ION XX(3),XXX(3) 

XX AND XXX ARE TOKEN DIMENSIONS WHICH WILL CHANGE AS THE 
DIMENSION STATEMENTS ARE CHANGED.. .THEY ARE USED FOR PICKUP RUN 
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EQUIVALENCE I VAR (1 1 • XVAR ) , (VAR( 2 1 »H) .( VAR( 31 ,C IC II I 
C VAR(2«-IMAX«-1)=EVI(1 ) 

EQUIVALENCE! M» XX (1 ) ) , { OELX, X XX ( 1 ) ) 

C 

REAL MUI,LAMSQ.MUINF,MU 
REAL N I 

REAL MINF.MINFSQ, LAMBDA 
REAL NUI J , NUPI J 
C 

C INPUT 

C 

NAME L IST/NAM1 / DELX , ZSTERM, IMAX, JMAX, IBUS ,MJ, M, R, GAMMA, 

1 CIINF,PINF»TINF,VINF,MUI» THETAI »DGENI , FI , DELH I, .DEL l. ,EVI ,LI, 

2 GIL*EPSIIL»AJ»BJ»EJ»AIJ»NUIJ»NUPIJ»ALPIK* BETA IK ,SI G IK, X I , Cl MAX, 

3 ele i,ele2,xpst,nxpst,ipf, timer, ipickup.hcheckt.tcheckt.phmax 
C 

EXTERNAL FOFTS 

C INITIALIZE 

REMIND S 
REMIND 9 
REMIND 10 
REMIND 11 
REMIND 12 

REMIND 13 
NERR=0 

C TIMER= ONE HOUR 

TI ME R= 3600.0 
DO 5 1=1,15 
5 I BUG ( I > =1 

NXPST=2 
I P IC KU P=0 
ICHEK=0 
ITNEG=0 
■I EXP-0 
NEG=0 
IGK=0 
IUNEG=0 

DELZS=. 001952125 
DELX=. 015625 
XI =• 000122070312 5 
Cl MAX=.0625 
IPF=5 

DO 10 1=1,52 
VARCI»=0. 

CUVARI I >=0. 

IF ( I.E0.52 ) SO TO 10 
DERI II =0. 

ELE1 1 1 1=0. 

ELE2 1 1 1=0. 

10 CONTINUE 

C READ INPUT 

READ (5,15) (INI I),T=1,8) 

15 FORMAT (8A10) 

READ ( 5 » NAM1 ) 

IF ( IPICKUP.EO.O I G3 TO 60 
PRINT 1$, (INI I ) ,1=1 ,8) 

C READ PICKUP FROM TAPE13 

C THIS SHOULD PICKUP ALL COMMON AND LABELED COMMON 

C PLEASE CHECK DIMENSIONS IF THERE ARE PICKUP PROBLEMS 

C 

READ (13) (XXI I ) ,1=1 ,8961) ,(XXX( I) ,1=1 ,128) 

DO 20 1 = 1,1 ZTE RM 

READ (13) X.ZS.RS.RC, COST, SINT 

CALL RECOUT (10, 1, 0 , X,ZS »R S, RC .COST, SI NT ) 

READ (13) (DPDX( J),J=I,IZTERM) 
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CALL R ECO UT (8,2,0,QP DX, I , IZTEfM.l I 
READ (13) IX, TS, ES . P S , RHOS ,US, PSIS ,HS 
READ (13) (EVIS(J) »J=1»IMAX) 

CALL RECOUT ( 1 1 , 1 , 0. I X , T S, ES , PS » RHOS ,U$ . PS I S , HS » 

CALL RECOUT ( 1 1 , 2, 0 , EV IS , 1 » I MAX, 1 1 
READ (13) X 

READ (13) (P(J ), J=I, IZTERM) 

CALL RECOUT 19,1, 0,X) 

CALL RECOUT (9 .2 ,0 . P , I , I ZTERM, 1 > 

20 CONTINUE 

IF (KSTAG.EO.O) GO TO 30 
DO 25 1 = 1 » KSTAG 

READ (13) K , RHOURT ,VAR(1 ) , PS IS 

CALL RECOUT (12,1,0, K, RHOURT ,VAR(1 ) » PSIS ) 

25 CONTINUE 

30 REWIND 3 

REWIND 9 
REWIND 10 
REWIND 11 
REWIND 13 
IREAD=IPSI 

IF ( IRSI.E0.6. AND, ISTAG.EO.l ) IREAD=IPSI-1 
DO 35 1=1. IREAD 

CALL RECIN ( 11 , 1 . 1 RE C I N, I X ,T S, ES , PS , RHOS ,U S , PS I S . HS ) 

CALL RECIN ( 11 ,2 , IREC IN, EV IS ,1 , I MAX, 1 ) 

35 CONTINUE 

DO 40 1=1. IREAD 

CALL RECIN ( 9, 1 , IR EC IN , X ) 

CALL RECIN (9,2. IREC IN, P, I, I ZTERM, 1) 

CALL RECIN ( 8. 2, IREC IN, DPDX, I , I ZTERM ,1 ) 

40 CONTINUE 

READ (131 (XX( II .1 = 1,8961) ,( XXX( II ,1=1.128) 

REWIND 13 
PRINT 45 

45 FORMAT ( 24HJ PICKUP CASE /> 

PRINT 50, IMAXL. (XPST(I ) ,1=1 ,101 .( ELE2 (I) . 1 = 1.5) 

50 FORMAT ( 1 X , 6HI MA XL = , 1 3 , / IX ,5HXPST= , 2 (5 E 17. 8/ » / IX , 5HE LE2= , 5 El 7, 8/ IX 
1 .6HKST AG=, I 3/ I 

PRINT 600, I PSI , IX ,X VAR ,H,MU,PFTL,RH0,U,T, E.SPEC , SUMCI 
IM AXP1 = IMA X+l 

PRINT 605. (VAR( I) ,t=3,IMAXPl) 

PRINT 610. (CM (I », 1 = 1, IMAX I 
IMAXP3=2+IMAX*1 

PRINT 615, (VAR( I) , i=IMAXP3, IMAXL) 

SPEC=0. 

IX=IX-1 
IBUG(1 1 ) = 0 
KEYINT =0 
K1PF=IPF-1 
HCHECK=HCHECKT 
TCHECK=TCHECKT 

PRINT 55, IPF.CIMAX, TIMER, HCHECKT, TCHECKT, IZTERM, KSTAG 
55 FORMAT (1CH IPF = I2,10H C IMAX=E15. 8 , 10H TI MER = E15. 8 , 10H 

1 HCHECKT=E15.8.10H TCHECKT=E15. 8/10H I Z TERM=I 2. 10H KSTAG=I3/ 

2 ) 

GO TO 560 
60 CONTINUE 

C 

IF ( IBUG(1 J.E3.1 I GO TO 65 
WRITE (6.NAM1) 

65 PRINT 70 

70 FORMAT (1H1) 

PRINT 15, ( IN ( I) .1=1,8) 

C 

C SUM OF Cl INF SHOULD BE 1.0 

SUM=0. 
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00 75 1=1, IMAX 

75 SUM=SUM+CI INF( I) . 

IF (SUM.EO.1.0) GO TO 85 

PRINT 80, SUM, ICIINFd 1,1 = 1, IMAX) 

80 FORMAT <34H SUM OF CIINF SHOULD BE 1.0. IT IS, £17.8, 8H CIINF=/5<5 
l El 7. 81) 

STOP 1 

C STOP 1 

85 CONTINUE 

HCHECK=HCHECKT 

TCHECK=TCHECKT 

C READ IN SPECIES AND REACTIONS 

READ 15,90) ( S PECI E ( I ) , I =1 , 1 MA X ) 

90 FORMAT (8A10) 

PRINT 95 

95 FORMAT I//) 

PRINT 100 

100 FORMAT (12H SPECIES/) 

00 110 1*1. IMAX 

105 FORMAT CI5,A10) 

110 PRINT 105, I , SPECI E( I ) 

PRINT 95” 

JMAX2=2*JMAX 

READ 15,90) (REACTU) ,1 = 1, JMAX2) 

PRINT 115 

115 FORMAT <14H REACTIONS/) 

J=0 

DO 120 1=1 » JMAX2 ,2 
J=J+1 

120 PRINT 125, J, REACT (I ),REACT( 1*1 ) 

125 FORMAT (I5,2A10) 

PRINT 95 

C PRINT INPUT 

PRINT 130, IPF 

130 FORMAT I17H PRINT FREQUENCY 3 13 ) 

PRINT 135, HCHECKT .TCHECKT »PHMAX 

135 FORMAT (10H HCHECKT=E15.8 ,10H TCHECKT=E15. 8, 10H PHMAX=E15.8) 
NXPSTP1=NXPST+1 

PRINT 140, XI,CIMAX*EL£1,ELE2,NXPST, (XPST( I) ,1=1, VXPSTP1 ) 

140 FORMAT C4H XI=E15.8/7H C IMAX=E 15. 8/6H ELE1 =2E17. 8/7C 7E17. 8/1 /6H El 
1E2=2E17.8/7(7E17.8/I/7H NXPST=I3/6H XPST=/15(7E1 7. 8/ ) /) 

PRINT 145, DELX,ZSTERM, TIMER, IP1CKUP 

145 FORMAT (6H DEL X= El 5. 8/ 8H ZSTERM=E15. 8/7H T IMER*E 15. 8/9H IPICKUP=I1 
1//) 

PRINT 150, IMAX, JMAX ,M,R, GAMMA, P INF, TI NF , VINF , I BUG 
150 FORMAT <6H IMAX=I2,8H JMAX=I2,5H M=U/3H R=E15.8/7H GAMM A=E1 5. 
18/6H P I NF= E15. 8/6H T INF=E15. 8/6H VI NF=E15. 8/6H I BJG= 11,1415/ ) 

PRINT 155 

15 5 FORMAT I5X.3HMUI12X, 6HTHET A I 9X ,5HDGENI 10 X, 2HFI 13X, 5H0ELHI 10X , 4HDEL 
1I11X,5HCIINF10X,3HEVI12X,2HLI/) 

DO 165 1 = 1, IMAX 
160 FORMAT ( 7E15.5 , D15 .6 , 1 5 ) 

165 PRINT 160, MUI ( I ), THETA I ( I ) , DGEN I( I ) ,F I ( II ,0ELHI(I), DEL I (I ) * Cl INF ( 
1 I ) «E VI ( I ) , LI ( I ) 

PRINT 170 

170 FORMAT { /6X , 2HMJ6X ,2 HA J16X ,2HB J16X , 2HE J/ I 
DO 180 J=1 » JMAX 
175 FORMAT 118, 3E18. 8) 

180 PRINT 175, MJIJ),AJ(J),BJ( J),EJ(J) 

PRINT 95 
IMAXP1=IMAX+1 
DO 225 1=1, IMAX 

PRINT 185. I,(GIL(L,I),L=1,IMAXP1> 

185 FORMAT (4H 1=15, 23H (GIL (L , I ) ,1=1, IMAX + 1 ) / 1 3 C 7E18. 8/1 I ) 
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PRINT 

190 « 

190 

FORMAT 

<4H 


PRINT 

195, 

19 5 

FORMAT 

<4H 


PRINT 

200, 

200 

FORMAT 

<4H 


PR INT 

205, 

20 5 

FORMAT 

(4H 


PRINT 

210, 

210 

FORMAT 

<4H 


PRINT 

215, 

215 

FORMAT 

<4H 


PRINT 

220, 

22 0 

FORMAT 

<4H 


PRINT 

95 

225 

CONTINUE 


I M= I MA 

X + l 


PRINT 

215, 


PRINT 

230 

230 

FORMAT 

( 2X, 


PRINT 

95 


I,(EPSIIL(L,I),L=1, IMAXPl) 

I=I5,26H (EPSIIUL, I) ,L=l,IMAX*l)/3(7E18.8/>) 
I » ( A I J( I,J) ,J=1 ,JMAX) 

1 = 15, 21H ( A I J ( I ,J), J=l, JMAX 1/5 (7E 18.3/1 I 

I • ( ALPIKIK. I) , K=1 , 1 MAX I 
1 = 15. 2 3H (ALP IK CK, I ) , K= 1 , IM AX ) /3( 7E 18. 8 / ) ) 

I, (BETA IK (K , I I » K=l, IMAX) 

1=15 ,2 AH (BETAIK(K.I) ,K=1 , 1 MA X I /3 ( 7 El 3. 8/1) 

I « ( S I GI K( K, I ) ,K=1 , 1 MAX) 

I =15.2 3 H ( S IG IK (K , I ) , K= 1, I MAX 1/3 ( 7E 18.8 /) I 
I « ( NU I J ( I , J I , J=i, JMAX I 
1 = 15 ,22H (NUIJ(JtI),J=l,JMAXI/5(7 El 8.8/11 
I « ( NUPl J( l » J I , J=1 , JMAX) 

1=15.2 3 H (NUPIJ(J,I >,J=i,JMAX)/5(7E13.8/)) 


IM, ( NUI J ( I M , J ) , J= 1 , JMAX l 
9HEND INPUT// ) 


C 

C DEFINE SHOCK GEOMETRY 

C AT EACH X FIND Z S.RS ,RC, COST ,S INT AND' STORE ON TAPE 10 

C 


CALL SHOCK G (DEL X, ZSTERM , I ZTERM ) 


C 

C IZTERMsNUMBER OF DELTA X INCREMENTS AS GENERATED IN SHOCKG 

IF ( NXPST-l.GT.I ZTERM) PRINT 235 

235 FORMAT (84H NXPST-1. GT. I ZT ERM, IT MUST BE .LE. FOR PHYSICAL SPACE C 
1 ALCULAT10NS - SEE DO 700 LOOP/) 

IF ( NXPST-l.GT.I ZTERM) STOP 301 

PRINT 415. I ZTERM 

IF ( IZTERM.GT.5D0I PRINT 240 

240 FORMAT (1X.50HIZTERM.GT.500 CHANGE DIMENSION OF P.OPDX AND VARI ) 
IF ( IZTERM.GT.500 ) STOP 301 


C STOP 301 

C FREESTREAM QUANTITIES 


SUM=0 

DO 245 1 = 1 , IMAX 
245 SUM*SUM+CI INF( I)/MUI (I) 

MUINF=1. /SUM 

RHO I NF=MUI NF*P INF/ (R*TINF)- 
AINF=S0RT (GAMMA*R*TINF/MUINF ) 

MINF =V INF/ A INF 

C FOR EACH I, SPECIE 

EINF=0 

DO 275 1=1 . IMA X 

. TEM=EXP(THETAI ( I l/TINF ) 

EVIINF(I) = (R*THETAI( m/(MUI(I >*(TEM-1.) )*FI (I I 
C FOR EACH L, REACTION LEVEL 

GSUM=0 
GESUM=0 
Lii=Lim 

IF (LI I.LE.20) GO TO 255 
PRINT 250. LI I 

,250 FORMAT ( IX *4HL II =, 13 ,2X, 28HA LEVEL IN LI ARRAY IS GT 20/2X.43HNEED 
1 TO CHANGE DIMENSION OF EPSIIL AND GIL » 

255 CONTINUE 

00 2 65 L=1 « L 1 1 
TEM1=EXP(-EPSIIL(L,I l/TINF) 

GSUM=GSUM*GIL(L, I)*TEM1 
GESUM=GESUM+GIL(L,I I *E PS 1 1 L ( L , I ) *TEM 1 
IF ( IBUG(1 I.NE.O ) GO TO 265 

PRINT 260. I,L.GSUM,GESUM.TEM1,EPSIIL(L.I).TINF.GIL(L,I) 

260 FORMAT (36H I , L . GS UM , GESUM ,TEM1 , EPSI IL ,T INF, GI L/ 2 15 . 6E18, 8/ > 
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26 5 CONTINUE 

EEIINF=R/MUI (I>*IG£SUM/GSUM) 

EIINF=l.5*R*TINF/MUI UI+FI (I )*R*TI NF/MUI ( I )+EVIINF< I H-EEI INF+OELHI 

im/Muim 

EINF=EINF+EIINF*CI INFII) 

IF ( IBUG<2).NE.0 I GO TO 275 

PRINT 270, I.EINF.TEM.THETAI (I ),TINF,MUI(I I.EVIINFI I ) ,GSUM,GESUM ,E 
1 El INF, EIINF.FI (I ), DELHI (I) 

270 FORMAT (34H I , El NF , TEM, THETA I , TI NF ,MUI ,E VI I NF/ 15 . SE1 8.8/ 33H GSUM.G 
1ESUM ,EEIINF,EI INF, FI.0ELHI/6EI8.8/) 

275 CONTINUE 

PRINT 280, MUINF.RHOINF, AINF,MINF,EI INF, El NF. EE I INF 
280 FORMAT (31H ♦ FREESTREAM QUANTITIES ♦♦♦/8H MU1NF*E15.8,10H R 

l H0INF=E15. 8,8H A INF=El 5. 8, 8H M INF=E15. 8. 9H eiINF*El5.8/7H E 

2 INF= El 5.8/9H EE 1 1 NF=/ ( 8E18. 8/ )) 

PRINT 285, (EVIINFII ), 1=1, IMAX) 

285 FORMAT < 9H EVI INF*/ <6E18. 81 I 

C END FREESTREAM QUANTITIES 

C BEGIN QUANTITIES BEHIND SHOCK 

C FOR RANGE _OF_X 

PRINT 290 

290 FORMAT I//33H QUANTITIES BEHIND SHOCK ♦♦♦/6X.2HTS15X,2HES16X, 

1 2HPS 16X ,4HRH0S 14 X, 2HUS 16 X, 4HPS I S 14X, 2HHS I 
IE0F=1 0 

IF (NXPST.E0.2) PRINT 310 
IF <NXPST.EQ.2I STOP 2 
NXPSTMl=NXPST-l 
ITK=0 

DO 385 I X= 1 , 1 Z TERM 
C 

C USE RECOUT TO WRITE 

C USE RECIN TO READ 

C 

CALL RECIN < 13 ,1 , 1 RECIN, X, ZS ,RS , RC .COST, SI NT ) 

IF (ENDFILE 10) 295,305 
295 PRINT 300. IEOF 

300 FORMAT <14H IN MAIN EOF 13) 

STOP 6 

C STOP 6 

305 CONTINUE 

310 FORMAT (78H NXPST=2 - SURELY SOME PHYSICAL SPACE VALUES ARE OESIRE 
10 ~ MAKE IT 3 AT LEAST) 

DO 315 IT = 2 .NXPSTMl 

IF < ABS(XPST(ITI-X).GT.O.OOOOOl) GO TO 315 
ITK* IT K+l 
315 CONTINUE 

TEM=SINT**2 

LAMBDA=RHOINF*VlNF*SINT 

OMEGA=PINF+RHOINF*VI NF**2*TEM 

DELTA* EINF+PINF/ RHOI NF+( VI NF**2*TEM) /2. 

LA MS 0 s LAMBDA** 2 
0MEGS0=0MEGA**2 
MINFSQ = MINF**2 
TEM=TE M*MI NFSQ 

TG = (TINF*(2.*GAMMA*T EM- (GAMMA- 1.0) )* < < GAMMA-l. 0 > *TEM+2. ) )/<<GAMMA+ 
11. )**2*TEM) 

IF ( IBUG(5 ).NE.O) GO TO 325 

PRINT 320, X, SINT, LAMBDA, OMEGA, DELTA, TEM.TG 
320 FORMAT <39H X, SINT, LAMBDA, OMEGA, DELTA, TEM, T3/7E17.8/) 

C 

C C5.002 LF-ITR1 

325 CONTINUE 

DELTTS=10. 
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E 1 = • 00 0 1 
E2=» 001 
MAXI =50 

CALL ITR1 (TG»0ELTTS»F0FTS,E1»E2»MAXI» I CODE I 
IF (ICODE.EQ.il PRINT 330 

330 FORMAT (34H MAXIMUM ITERATION EXCEEDED *****//) 

IF (ICOOE.EO.2I PRINT 335 

335 FORMAT I30H DERIVATIVE =0. IN ITR1 ♦****//) 

IF ( ICODE.EQ.l.OR.ICODE.EQ.2) STOP 321 
C STOP 321 

TS=TG 
ES=0. 

DO 360 1*1 , IMAX 

IF (M.EO.l) EVISU ) = EVIINF(I ) 

IF ( M» EQ« 1 I GO TO 345 
TEM=EXP(THETA1(I )/TS> 

EVISII )=(R*THETAII I » » / ( MUI (I)* IT EM-1.) )*Fl (I ) 

IF ( ieUGm.NE.OI GO TO 345 
PRINT 340, IX, I » TEM , EV IS IT ) 

340 FORMAT (19H I X , I , TEM , EV IS 1 1 ) = 2I 5, 2E 18.8 ) 

345 SUMG=0. 

SUMGE=0. 

LI I *LI 1 1 • 

DO 350 L=1 ,L 1 1 

TEM1 =E XP (-EPSI IL IL » I )/TS> 

SUMG=SUMG+TEM1*GILIL,I ) 

350 SUMGE* SUMGE+TE Ml *G I L (L, I )*EPSI IL(L,I ) 

EE IS=I R/MUI 1 1 ) )MSUMGE/SUMG) 

El S=(1.5*R*TS)/MUI ( I IMF II I > *R*TS I /MUII I I ♦ 1.0*EV IS! I l+EEI S+DELHI ( I 
1 » / MU I C I ) 

IF I IBUGI5 I.NE.O ) GO TO 360 
PRINT 355, 'EEIS.EIS, SUMS, SUMGE 
355 FORMAT (22H E El S, E I S, SUMG ,SUMGE=4E1 8. 8 ) 

360 ES = ES«-EIS*CIINF( I) 

PS*SQR T ( OMEGSQ-2. *LAMSQ* (DELta-ES I ) 

RHOS*!MUINF*PS)/IR*TS) 

US=VlNF*COST 

PS I S=I RHOI NF*VINF*RS**2 ) /2. 

HS=PS/RHOS+ES 

C 

C WRITE ON TAPE 11 AT EACH X 

CALL RECOUT (1 1 , 1, 0, IX,TS, ES ,PS,RHOS ,US,PS IS ,HS ) 

CALL RE COUT 1 1 1 , 2. 0, EV 1 S . 1 , 1 MAX, 1 ) 

C ’ ” ’ . " ■ 

PRINT 365, TS»ES»PS , RHOS »US,PSIS,HS 
365 FORMAT I/7E18.8I 

PRINT 370, IX, (EVISII), 1*1, IMAX) 

370 FORMAT I5H IX*I3,7H EVIS*/4I 7E18.8/I /) 

IF I IBUGI5 I.NE.O I GO TO 385 
PRINT 375 . 

375 FORMAT I50H /////////////////////////////////////////////////) 
PRINT 380. IX, TS ,P S.RHOS , US, PS IS«HS,ES 
380 FORMAT I29H I X, TS ,P S,RHOS ,US, PS IS ,HS, ES/ 1 4,7E 18. 8/ I 
385 CONTINUE 
REWIND 10 
END FILE 11 
REWIND 11 

PRINT 390, ITK.NXPST 

390 FORMAT I / /2X, 1 OHXX XX X I TK= 13 ,2 X , 6HNXPST* 1 3 , 2X . 5H XXXX X// ) 

IF ( ITK.GE.3*NXPST/4) GO TO 400 
' PRINT 395 

395 FORMAT I 89H NUMBER OF X-S IN XPST ARRAY. I TK, IS .LT. 3*NXPST/4 - R 
1EEXAMINE DELX, NX PST AND XPST ARRAY) 

C SECOND STOP 301 

STOP 301 
400 CONTINUE 
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c 

i c 

c 

, c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 

_c 

c 

c 


c 


40 5 
410 
C 


c 

415 


C 

c 

420 

425 


C 


43 0 
43 5 


440 


END COMPUTATION OF QUANTITIES BEHIND SHOCK 


CALL FOPOX 

FIND PRESSURE, P, AT EACH X FOR EAC-I STREAMLINE AND 
STORE ON TAPE 9 

FIND P DOT, DPDX, AT EACH X FOR EACH STREAMLINE AND 
STORE ON TAPE 8 

INTEGRATE ALONG EACH STREAMLINE PS I 
VAR ( 1 ) = X VAR 

VAR<2>= H, ENTHALPY VARC2I MUST BE H, WHICH MAY BE + OR - 

CONCENTRATION OF SPECIE 

VA R( 3 ) = CI(1I 

VAR(2+IMAX)=CI( IMAX) 

EQUILIBRIUM VIBRATIONAL ENERGY 
VAR ( 2* I MAX + 1 ) = EV I ( II 

VAR(2+IMAX + IMAXI=EVI UMAX) 


DO 410 1=1, IMAX 
TRUNCATE 

Y=F I (I ) * ( DEL I ( I ) /THE TAIC I ) +1.0 ) 
NI(I)=AINT<Y> 

IF ( IBUGC5UNE.0 I GO TO 410 
PRINT 405, I , Y,N I( I ) 

FORMAT (12H I , Y ,N I (I) =1 5 , 2E IB. 8/) 
CONTINUE 


MU=0. 

N= 1 + 2* I MAX 
IF (M.EQ.O) N=N- IMAX 

INITIALIZE 

PRINT 415 * IZTERM 
FORMAT ( 9H IZTERM=I3) 

KSTAG=0 
ISTAG=0 
XP ST ( 1 ) = DE L X 

XP ST < N XP ST ) = VAR T ( I ZT ER M > +1 00*0 
SUMC 1=0. 

PRINT 95 
IZTM1= I ZTERM-1 
LPS1 =1 

IPSI=0 

3 EG IN EACH STREAMLINE COMPUTATION HERE 

IPSI=IPSI+1 

IF ( IPSI.EQ. IZTERM) GO TO 735 
CONTINUE 
K I PF =0 
I EOF =1 1 

ISTAG = 2 FOR STREAMLINE DELX 
IF (ISTAGoEQ.il 1ST AG=2 

CALL RECIN ( 11 ,1 . IRECIN, IX ,T S, ES ,PS ,RHQS ,US .PS IS.HS > 

CALL RECIN ( 11 ,2 , 1 RECI N, EV IS , 1 , I MAX, 1) 

IF (ENCFILE 11) 430,440 
PRINT 435, I EOF, IP S I 

FORMAT ( 5H EOF 1 4, 57H IN MAIN, TRY TO COMPUTE PHYSICAL SPACE VALU 

1 ES I PS 1= 13 // ) 

GO TO 735 
CONTINUE 
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c 

IF ( IPSI.EQ.l) H ST AG = HS 
C 

CJ=XI 

SPEC=0.0 

C EVALUATE DERIVATIVES WHEN SPEC=0.0 

11=0 

C USE SHOCK VALUES FOR INITIAL COMPUTATION ON EACH STREAMLINE 

KEYINT=0 
IX=0. 

T=TS 

IF (M.EQ.O) GO TO 450 
DO 445 1=1* 1 MA X 
445 TV I C I) =T I NF 
GO TO 450 

450 DO 455 1*1, IMAX 
'455 TVim-TS 
460 RHO= RHOS 
MU=MUI NF 
U=US 

DO 465 1=1, IMA X 
MM*2+I 
K= IMAX+MM 
VARIMM)=0. 

VAR(MMI=CIINF( II 
VAR I K) =0. 

IF (M.EO.Ot VAR( K) =E VI S( It 
IF (M.EO.l) VAR(K)*EVIINF(I) 

EVII I ) =VAR ( K I 
465 CONTINUE 
E=ES 
H=HS 

IF I IPSI.E0.6, AND, ISTAG.EQ.OI GO TO 635 
C 

I E OF =9 

CALL RECIN ( 9, 1 , IREC IN , X ) 

CALL RECIN (9, 2 , IR EC IN, P , I PS I , I ZTERM ,1 1 
IF IENDFILE 91 430,470 
C 

470 IE0F=8 

CALL RECIN (8, 2, IR EC IN,DPDX, IP SI ,1 ZTERM, 1 1 
IF (ENOFILE 81 430,475 
475 IF (JBUG(3).NE-0) GO TO 485 

PRINT 480, (VARI (IT) ,DPDX( IT ) , P( IT >, IT=I PS I . IZTE RW » 

480 FORMAT (33H ( DPDX ( I T ) , P ( I T) , I T= IPSI , I ZT ERM) / ( 6E 17. 8 I) 

485 CONTINUE 
VARI 1 ) =X 

C OMIT PSI=0.0 STREAMLINE FOR NOW PICK IT UP LATER 

IF (X.NE.O. ) GO TO 490 
IG=0 

PRINT 725 
GO TO 730 
C 

490 CONTINUE 

IF ( X.EQ.VARI ( IPSI ) I GOTO 500 
C 

PRINT 495, IPSI, X, VAR H IPSI) 

495 FORMAT <6H IPSI=I4,4H X=E15.8,13H VARI 1 1 PSI ) =E15. 0 ,34H X AND VA 
lRI(IPSI) SHOULD BE EQUAL/63H EXAMINE GENERATION OF X,ON TAPE 9, AN 
2DVARI IN FDPDX SUBROUTINE) 

STOP 663 

C STOP 663 

500 CONTINUE 
C INITIALIZE 

C 
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EL B= 0 
SPEC=0 

IF ( VAR(1).GT.XPST(LPS1)> LPS1*LPS1+1 
LP S2=L PS1 
505 CONTINUE 

C CHECK ON PICKUP 

CALL SECOND (CHECKT) 

IF ( IX •NE, 1 ) GO TO 515 
PRINT 510 * CHECK T 

510 FORMAT (20H -CHECKT =E15# 8 ) 

515 IF (ICHEK. E0*1) GO TO 540 

IF ( ( T IMER-CHECKT) #GT# 30# ) GO TO 560 
IF < IP SI«GT#5 ) GO TO 525 
PRINT 520 

520 FORMAT ( IX • /10X, 1DHSS$$S$$$$S/67H MUST BE ON OR BEYOND PSI=0«, 

1 I PS 1 = 6* BEFORE PICKUP IS POSS IBLE/10X f 10H$$$$ $$ %%%% / ) 

525 CONTINUE 
K I PF =0 
I P F= 1 

C WRITE PICKUP TAPE 

WRITE (13) (XX(n»I=l f 8961),(XXXJI)*J = lf 128) 

PRINT 530, TIMER, CHECKT 

530 FORMAT (30H PICKUP TAPE WRITTEN — T IMER=E 15« 8 , 1 OH CHECKT= E15. 8/ 

1 ) 

I CHE K= 1 

PRINT 600, IPSI, IX,XVAR,H,MU,PFTL,RHO,U,T, E,SPEC,SUMCl 
PRINT 605, (VAR( I) , I =3 , 1 MA XP2 ) 

PRINT 610, . (CM ( I ) f 1=1, IMAX) 

IMAXP3=2+IMAX+1 

PRINT 615, (VAR( I) ,I=IMAXP3, IMA XL) 

PRINT 620 
PRINT 535 

535 FORMAT (/93H — 

ITEST=0 

540 IF (TIMER-CHECKT.GT. 10*) GO TO 560 
REWINO 8 
REWIND 9 
REWIND 10 
REWINO 11 
REWIND 12 
DO 545 1=1,1 ZTERM 

C READ 

C WRITE ON 13 TO SAVE FOR PICKUP 

CALL RECIN ( 10 ,1 , I RECI N , X, ZS , R S, RC ,C OST, SI NT ) 

WRITE (13) X,ZS,RS,RC, COST, SINT 

CALL RECIN ( 8, 2, I R EC I N , DPDX, I , I Z TE RM . 1 ) 

WRITE (13) (DPDX(J),J=I, IZTERM ) 

CALL RECIN ( 11 , 1 , 1 REC I N, I X ,T St ES , PS , RHOS ,U S , PS I S , HS ) 

CALL RECIN ( 11 , 2 , I REC IN , E V IS , 1 1 1 MAX, 1 > 

WRITE (13) IX,TS,ESfPSfRHOS,US,PSIS,HS 
WRITE (13) ( EV IS ( J ) ♦ J=1 , IMAX) 

CALL RECIN ( 9, 1 , IR EC IN , X ) 

CALL RECIN (9, 2, IREC IN , P , I , I ZTERM* 1 ) 

WRITE (13) X 

WRITE (13) (P( J) ,J=I , IZTERM) 

545 CONTINUE 

DO 550 1=1 , KST AG 

CALL RECIN ( 12 , 1 , I REC I N, I P SI , R HO URT , VAR( 1 ) , P SI S > 

WRITE (13) IPSItRHOURT ?VAR(1 ) , PS IS 
550 CONTINUE 

PRINT 50. IMAXL, (X PST (II ,1=1,10) , ( EL E2 ( I > , 1 = 1, 5) ♦ KST AG 
PRINT 555 

555 FORMAT (25H DATA ON TAPE 13 /) 
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STOP 422 

C PICKUP DATA ON TAPE 13 STOP 422 

560 CONTINUE 
IX=IX+1 
C 

C CALINTH, NOD IF I E 0 RUNGE-KUTTA 

CALL CALINTH ( VA R, DER t EL El , ELE 2» CJ . S PEC, N , CU VAR , ELB , Cl MAX , NE RR ,PHM 

1 AX ) 

C 

IF ( NERR.EO.O) GO TO 575 
IF (NERR.EO.l) PRINT 565, CJ,N 
565 FORMAT ( 16H BAD INPUT, CJ ,N=E 17.8 , 1 5 ) 

IF (NERR.E0.2) PRINT 570, ELE1 ,ELE2 
570 FORMAT ( 1 X ,21HEX AMINE ELE1 AND ELE2/1GI7E17.8/ ) > 

IF (NERR.EO.l. OR. NERR.EQ. 2 ) STOP 665 

C STOP 665 

575 CONTINUE 

IF ( IX. EO. 1> GO TO 580 

IF ( IPSI.EQ.6. AND. IX.EQ.6) GO TO 580 

IF ( VAR ( 1 ) • GE» VARI ( I ZTERM) I GO TO 580 

IF ( IPSI.LT.6.AND.VAR(1).GE.DELX» GO TO 580 

KIPF=KIPF+1 

IF (KIPF.NE.IPF) GO TO 625 
KIPF=0 
C 

580 CONTINUE 

IM AX P2 = I MAX + 2 
SUMC 1=0. 

DO 590 I SUM=3, IMAXP2 
IF (VAR( ISUMI.LT. 0. I NEG=-1 
IF (VAR(ISUM).LT.O. > PRINT 595 

585 FORMAT (/5X.21H NEGATIVE Cl- /) 

590 SUMC 1= SUMC I ■♦•VAR ( I SUM ) 

DO 595 ICM=1 , 1 MA X 
C 

C PRINT ANSWERS 

595 CM(ICM) = VAP(ICM+2)*MU/MUK ICM) 

PRINT 600, IPSI, IX,XVAR,H,MU,PFTL,RHC,U,T,E,SPEC,SUMCI 
600 FORMAT (2X12,15. 8X ,2 HX=D 15. 8 ,7X , 2HH=D15. 8 , 6X , 3 HMU=E1 5. 8, 7X , 2HP=E 15 
1.8,5X,4HRH0=E15. 8/17 X ,2HU=E1 5. 8 , 7X ,2 HT =E1 5 . 8 , 7 X . 2HE= El 5. 8 , 6X , 3HC 1 = 

2 El 5. 8, IX, 8H SUM (CIt=E15«8) 

PRINT 605, (VAR( I ) » I =3 » I MA XP2I 

605 FORMAT (10X.9H C I ( 1)=D15.8,9H C I ( 2)=D15.8.9H CI( 3)=D15.8,9H 

1 CI< 4 ) =015.8. 9H C I ( 5 ) =015. 8/10X *9H CI( 61=015. 8. 9H C I ( 71=015 
2.8.9H CI( 8 ) = Dl 5, 8, 9H CI( 9)=D15.8,9H C I ( 10 ) = D15. 8/lOX, 9H CI-(1 

311=015. 8, 9H CK 12» = D15.8,9H Cl ( 13) =D15.8 ,9H C I< 14 I =01 5. 8 , 9H ' Cl 
4(15)=D15.8/1CX,9H C I ( 16 )=D15. 8, 9H Cl ( 17) =015. 8. 9M Cl ( 18 ) =015. 8, 
59H Cl <19)=D15.8,9H C I ( 20 ) = D1 5. 8/ 10 X, 9H Cl ( 2 1) =015. 8. 9H CI(22)= 
6015.8, 9H Cl (23 ) =015.8, 9H C I( 24 )=01 5. 8, 9H Cl (25) =D15.8) 

PRINT 610. ( CM( I ), 1=1 , IM AX ) 

610 FORMAT (10X.9H CM ( 1>=£15.8,9H CM( 2)=E15.8,9H CM( 3)=E15.8,9H 

1 CM( 4 )=E15,8,9H CM( 5 ) =E 15. 8/lOX ,9H CM( 6)=E15.8,9H CM( 7)=E15 
2.8.9H CM( 8 ) = E15. 8..9H CM( 9I=E15.8,9H CM( 10>=E15. 8/10X.9H CM(1 

31)=E15.8,9H CM( 12 )=E15. 8,9H CM( 13) =E15.8 ,9H CM( l 4 )=E15. 8, 9H CM 
4(151 =E15.8/10X,9H CM( 16 ) = E15. 8, 9H CM( 1 7 ) -E15.8 . JH CM( 18 ) = E15. 8, 
59H CM(19)=E15.8,9H CM{ 20 )=E15. 8/ 10 X, 9H CM(21I =E15.8,9H CM (22) = 
6E15.8.9H CM( 23 ) =E15.6 ,9H CM( 24 ) =E15. 8, 9H CM ( 25 1 =E 15. 8 ) 
IMAXL=2+IMAX-MMAX 
IMAXP3=2+IMAX+1 

PRINT 615, (VAR( I) ,I=IMAXP3, IMAXL) 

615 FORMAT (10X.9H EVI( 1)=015.8,9H EVI( 2)=D15.8,9H EVI ( 3)=015.8,9H 
1 EV I ( 4 ) =015.8, 9H EVK 5 )=015.8/10X,9H EVI( 6)=D15.8,9H EVI( 7I=D15 
2.8.9H EVK 81=015.8, 9H EVI( 91=015. 8. 9H EV I ( 10 » =015. 8/10X.9H EVI (1 
311=015. 8,9H EVK 12 1=015.8, 9H E VI ( 1 3) =015. 8 ,9H EV I ( 14 ) =015. 8, 9H EVI 
4(151=015. 8/10X.9H EV I ( 16 ) = 015. 8, 9H EVH17.I =015.8, 9H EVK 181 = 015. 8, 
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S9H EVI(l9)=D15.8,9H EVK20 1 = 015. 0/ 10 X* 9H E V I( 21 ) = Di 5. 8* 9H EVI<22)= 
6D15.8, 9H EVI(23I=D15.8,9H EV l< 24 ) = D1 5. 8, 9H E VI I 2 5 1 =015. 8 ) 

PRINT 620 
620 FORMAT {/) 

C 

IF ( IBUG(11 I.EO.O) I BUG< 11 ) = 1 
625 CONTINUE 

IF (NEG.EO.-l) PRINT 630 

630 FORMAT (31H A NEGATIVE CONCENTRAT ION, Cl I 

IF (NEG.EO.-l) STOP 670 

C STOP 670 

C 

C WHEN I P$ I =6 » DEL X AND STAGNATION STREAMLINES WILL BE COMPUTED 

C 

C I UN EG = 1 IF U**2 NEG, IN BASIC 

IF ( IUNEG.EO.l > GO TO 735 

IF ( IPSI.LT.6.AND.VARm.LT.XPST(in GO TO 505 
IF ( IS TAG. EO.l.OR. I STAG. EQ.2 ) GO TO 690 
C I ST AG=0 UNTIL AFTER EXtRAPOLATI ON FOR P$I=0. STREAMLINE 

C I ST AG= 1 FOR PSI =0.0 STREAMLINE 

C ISTAG=2 FOR DEL X STREAMLINE AND THEREAFTER 

C 

635 CONTINUE 
I G= I G+l 
IGC=0 

IG8= IMAX+2 
PS IG ( I G ) = PS l S 
DO 640 IG1=3,IMAXP2 
IGC= IGC+1 

CIG( IGC * IG)=VAR< IGl ) 

IGE= IGE+1 

EIG( IGC * IG ) =VAR( IGE ) 

640 CONTINUE 
C 

EG ( I G) =E 
TG(IG)=T 
HG(IG)=H 

IF ( IPSI.LT.6) PRINT 725 
IF ( IPSI.LT.6) GO T3 730 
C I PSI=6 STAGNATION STREAMLINE 

C AT DELX EXTRAPOLATE FOR PSI *0* VALUES 

PSI=0. 

MG=1 

CALL FTLUP (PSI .-Et MG t IGtPSIGt EG) 

CALL FTLUP ( PS I , T, MG , IG, PS IG,TG) 

CALL FTLUP ( PS I ♦ H, MG ♦ I G, PS IG ,HG ) 

C 

DO 650 IG2 = 1 # I MAX 

DO 645 IGl =1 ♦ 6 

CIT< IGl) =C IGC IG2.IG1 ) 

645 EIT( IGl ) =E I G( IG2,IG1> 

CALL FTLUP (PSI « Cl (IG2 ) » HG t IG« PS IGtCIT ) 

IF <CI(IG2).LT.0.0) Cl (I G2 1=0.00000001 
CALL FTLUP (PS 1. EVH IG2) tMG. IGtPSIGt EITI 
13=2+1 MAX+ I G2 

PRINT 646 » IGtI3t CITtCI (IG2 ) t EITt EVII IG2)tPSIG 
646 FORM AT (2I5/C7E17.8M 
650 VAR( 2+ IMAX+ IG2 )=EV I ( IG2) 

DD 655 IGl =1,6 

CALL RECIN ( ID ,1 , 1 REC I Nt X, ZS t RS, RC .COST, SI NT ) 

655 CONTINUE 

C INITIALIZE FOR PSI=0. STREAMLINE 

CJ = XI 
SPEC=0.0 
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11=0 

KEY I NT = 0 
IX = 0 
SUM=0 

DO 660 1=1 t IMAX 
SUM= SUM+C I \ I>/MUI( I ) 

IF (EVKI J.E0.0. ) TV I ( I ) =0# 

IF (Evnn.EQ.OJ GO TO 660 

XL=( (OGENU I ) * R*TH ET AI ( I )) /( MU I ( I) *E VI ( I)) ♦ 0 I 

ALN=ALOG(XU 
TVI II) =THE T AI ( I) /AIM 
660 CONTINUE 

PRINT 661 • ( EVI (I ) « TV I ( I ) « I S 1 « I MAX I 
661 FORMAT (5X*D17«8*E17#8I 
C 

U=SORT ( 2# 0*lHSTAG-H) ) 

C 

VAR ( 1 ) =0E L X 

C COMPUTE P AND DPCX AT IPSI=6 

PRINT 665 

665 FORMAT (24H PSI=0.000 STREAMLINE — ) 

P ( 6 ) =PS + US/(RC*R$) M-PSIS) 

RHO= P ( 6 ) / C H— E ) 

DO 670 1=7 v IZTERM 

CALL RECIN ( 11 • 1 # I RECI N # I X ,T St ES tPS t RHOS • US* PS I S #HS ) 

CALL RECIN 111 t2 ♦ I RE C IN* EV IS *1 1 1 MAX. 1 ) 

CALL RECIN ( 10 • 1 • I RE C I N, X * IS # R S, RC •COST, S I NT ) 

670 P< I)=PS+US/(RC*RS) *<~PSI$) 

REWIND 10 
REWIND 11 
MG=1 

NPG= IZTERM- 5 
DO 675 I =7 » I ZTER M 
LG= I -6 

675 DPDX(I I =01 F ( LG #MG* NPGt VARI (6)*P(6) I 
DO 6 80 1=1 1 5 

CALL RECIN (11 «1 * l RECIN# IX # T $# ES *P S • RHQS #U $#PS I $ # HS ) 

CALL RECIN ( 11 #2 # I REC IN# EV IS # 1 # I MAX* 1) 

680 CONTINUE 
I $ T A G= 1 
PSIS=0. 

CALL SECOND (CHECKT ) 

P RINT 6 95 * CH ECK T 

685 FORMAT (14H -- CH ECKT =E 15*8 #6H /) 

PRINT 686* (EVI (I lf 1 = 1# IMAX ) 

686 F0RMATI7D17.8) 

GO TO 500 
C 

690 CONTINUE 
C 

C AT EACH X WHERE PHYSICAL SPACE CALCULATIONS ARE DESIRED SAVE 

C 1. /(RHO*U> ON EACH STREAMLINE 

C 

IF ( VAR(1).LT.XPST(LP$2) > GO TO 720 
IF ( VAR(1).GT.XPST(LPS2> ) GO TO 705 
C V AR ( 1 ) = XPST1LPS2) 

RH0URT=1./(RHD*UI 
IF ( IBUG(13).NE#0) GO TO 700 
PRINT 695* I PS I * RHOURT 
695 FORMAT (7H IPSI=I3*9H RHOURT = E 15« 8 > 

700 CONTINUE 

CALL RECOUT ( 1 2* 1 * 0 * I PS I ,RHOUR T, VAR ( l ) * PS I S ) 

GO TO 715 

7C5 YPSI S= PREP SIS* (XPST( LPS2 )-PREX)*(( PS I S-PRE PS I S ) / ( VAR ( 1 ) -PREX H 

RHOURT = PRERU-M XPST ( L PS2 ) -PREX) ♦( ( R HO *U-PRE RU ) / (VAR( 1 >-PREX ) ) 
RHOURT = 1. /RHOURT 
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CALL RECUUT (1 2 * 1 . 0 . I P S I , R HOUR T , XP ST ( L PS2 ) , YPS I S ) 

IF ( I8UG(13).NE.0> GO TQ 715 
PRINT 695 ♦ I PS I • RHGURT 

PRINT 710 ♦ LPS1,LPS2,VAR(1),XPST(LPS2I.PRERSIS,PSIS,PRERU,P9EX 
710 FORMAT (12H L PS 1, LP S2 = 2 I 4 , 3 3H X, XP ST ♦ PRE PS I S , P SI S . PRERU , PR EX =/ 01 
17. 8. 5E 17. 8) 

715 CONTINUE 

LPS2=LPS2+1 
KSTAG=KSTAG+1 
720 CONTINUE 

PR EX=V AR ( 1 ) 

PRFRU=RHO*U 
PR EP SI S = P S I S 

IF I VAR ( 1 I .LT. VAR I ( I ZTERM) ) GO TO 505 
C FNO OF A STREAMLINE 

PRINT 725 

125 FORMAT </20H X< XX XX XXXX XX XX XX XX / ) 

IF i lPSl.EG.6o ANO. ISTAG.EG*1 ) GO TO 425 
730 CONTINUE 
GO TO 420 
C 

C COMPUTE PHYSICAL SPACE VALUES 

C 

735 CONTINUE 

CALL SECOND (CHECKT) 

PR INT 635 ♦ CHECK T 

REWIND 8 

REWIND 9 

REWIND 10 

REWIND 11 

PRINT 740. K ST AG 

740 F0RMATI37H REGIN PHYSICAL SPACE CALCULATIONS AT.I5.10H LOCATIONS/) 
C 

DO 840 LP$=1.NXPST 
REWIND 12 
I K =0 

745 CALL RECIN ( 10 . 1 . 1 RE C I N» X. ZS .R S , RC .COST. $ I NT ) 

IF (ENOFILE 10) 845.750 
750 CONTINUE 

IF ( IBUGI 13J.NE.0) GO TO 760 
PRINT 755 . LPS.X.XPST(LPS) 

755 FORMAT { 18H L PS . X , XPST ( LP S ) I5.2E30.15/) 

760 CONTINUE 

IF ( ABS(X-XPST<LPSn.GT.1.0E-8 ) GO TO 745 
IF ( IBUGI 13).NE.O) GO TQ 765 
. PRINT 755. LPS . X . XPST ( LP S ) 

C 

76 5 CONTINUE 

00 780 IPS = 1.KSTAG 

CALL RECIN ( 12 , 1 , I RE C I N , I P S I , RHO URT , VAR (1) . PS I S ) 

VARMX=VAR(1 )-XP$T(LPS> 

IF < ABS(VARMX).GT.l«0E-8> GO TO 780 
IK= I K+ 1 

IF (IK.GT.200) PRINT 770, IK 

770 FORMAT 1 1 X .50HIK.GT. 200. CHANGE DIMENSION OF RUT AND PSISTG, IK = .1 
13) 

RUTI IK ) -RHOURT 
PS I STG ( I K ) =PSI S 
IF ( IBUGU3UNE.0) GO TO 780 
PRINT 775, PSIS.RHOURT 
775. FORMAT (14H PSI S, RH0URT=2 E17. 8) 

780 CONTINUE 

C FIND SMALLEST DELTA PS I 
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$MALL=200o 
00 785 I =2 ♦ I K 

IF (PSISTGI II-PSlSTG(I-i)«LT •SMALL! SM ALL* PSISTGUI - PSISTGI 1-1 ! 

I 785 CONTINUE 
] XDEL=PS I STG< IK ) /SMALL 

C TRUNCATE 

L= I NT( XDEL+1# ) 

C REMAINDERING M*0 FOR EVEN* M=1 FOR ODD 

M=MOO< L.2 ) 

C MAKE L EVEN 

IF (M.NE.O! L = L + 1 

C FIND INTEGRAL 1«/(RH0*U) DELTA PSI FROM BODY TO SHOCK 

C SIMPSON-S RULE L INCREMENTS* L+l POINTS 

FL=L 

DP S I S s PSI STGCI K! /FL 

psii=o* 

RB-RUT < 1 ) +RUTI IK) 

LM 1=L- 1 

IF ( IBUG(i3I.NE«0) GO TO 795 
PRINT 790* L* M* SMALL *XDEL * DPSI S«R8 
790 FORMAT <25H L * M, SMALL * XDEL ,DPSI S*RB/2 15 , 4E17. 8 ) 

795 CONTINUE 

DO 800 I = 2 * L * 2 
PS 1 1 =PS 1 I +-DPSI S 

call ftlup (psii *R ui,it ik,psistg*ruti 

IF (I.EO.U RB = R8+4**RU1 
IF ( UEQ.L) GO TO 800 
PS 1 1 =P S 1 1 + DPSI S 

CALL FTLUP ( PSI I *RU2 ,1* I K, PS ISTG*RUT ) 

RB=R B+4« *RU1>2« *RU2 
800 CONTINUE 
C 

RB=RB«DPSIS/3. 

IF < IBUGU3).NE.0) GO TO 810 
PRINT 8C5 * RB*RS*COST*RB*$INT 
805 FORMAT (21H R B* RS ♦ COST * RB , S INT=5E17.8 I 
810 CONTINUE 

ARR= SORT IR S** 2-2 •* COST *R B) 

YI = (RS-ARR)/C0ST 
ZE=ZS+YI*$INT 

PRINT 815* X*ARR*Yl*ZEtPSISTG(l) 

815 FORMAT < 3H X=E15.8 *5H R=E15 # 8,5H Y=E15.8*5H Z = E15.8,7H PSI 

JL=E19.8) 

C " 

DO 830 I =2 • IK* 1 
DPSI=PSISTG(I)-PSISTG( 1-1! 

C TRAPEZOIDAL RULE 

TR = ( OP S I/2« )*(RUT< I ) +RUT <1—111 
R8 = R B- TR 

IF ( IBUG(13I*NE.0) GO TO 325 
PRINT 820* TR * RB 
820 FORMAT (5HTR* RB* 2E 17* 8 I 

82 5 CONTINUE 

ARR-SQRT (RS**2-2.*COST*R3> 

YI =( RS-ARR ) /COST 
ZE = Z S + YI *S I NT 

PRINT 815* X * ARR * Y I » ZE * P S l STG( 1 1 
830 CONTINUE 
PRINT 835 

83 5 FORMAT </! 

C ■ 

840 CONTINUE 

845 CONTINUE ^ 

PRINT 850 

850 FORMAT <29H END THIS CASE* HALLELUJAH) 
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STOP 

end _ 

subroutine basic 

BASIC CALLED BY CALINTH TO EVALUATE DERIVATIVES H DOT, 

C SUB I DOT AND EV SUB I DOT 
DERIVATIVES START IN DER(l) 

COMMON M,IX,IMAX,IPSI,IBUG(15) 

COMMON P< 500 1 »0.PDX( 500t,VARI( 5001 

COMMON EVIINF(25>,THETAI(25),MUI(25),FI(25),DELHI(25),CIINF<25), 
1 L I (25) 

COMMON R.MUINF, DELTA, LAMSQ*OMEGSQ. OMEGA, MU, T,U 
COMMON EPS I IL< 20 ,251 » GIL 120* 25 ) 

COMMON VAR (52) «CUVAR(52I,0ER(51) 

COMMON ZS.RS,RC, SI NT , COST, PS , US, PS I S 
COMMON MJ ( 50 ) , E J (50) , AJ ( 50 ) , BJ (50) / 

COMMON TVH25) ,N K 25 » ,DGEN I( 25 ) 

COMMON NUIJ(26,50I.AIJ(25,50).NUPIJ(25.50) 

COMMON SIGIK(25,25),ALPIK( 25,25) ,BET AIM 25,25) 

COMMON /LAB2/ DELX.ZSTERM, 12 TERM 

COMMON /LAB3/ El NF , P I NF , RHOI NF , V INF , E. JM AX , KEY I NT . RHO.HSTAG 
COMMON /LAB4/ PF TL ♦ K ITR1 , E ITR1 
COMMON /LAB5/ IUNEG 
COMMON /LAB7/ ITNEG.IEXP 

DOUBLE PRECISION VAR 

DIMENSION PHI ( 50 I , KJ (50) , S J( 50 I 

DIMENSION Cl (25) ,EVI (25) ,DCIDX(25) ,OEVIDX( 25) ,EVIBAR(25) 

EQUIVALENCE (CUV AR ( 2 ) , H) , ( CUVARI 3 ) ,C I ( 1 ) ) 

EQUIVALENCE ( DER (1 ) , DHDX ) , (DER (2 ) , DC IDX( I) ) 

REAL K J, LA MSQ, MU ,MUI , MU I NF ,N I, NU I J ,NUP I J 

DER ( 1 ) =DHDX MUST BE DHDX AS H MAY BE * OR - 

DER (2 )= DCIDX(l) 

DERU+IMAX l=DCIDX( IMAX) 

DER ( 1 + 1 M AX+1 ) *DE V I DX ( 1 ) 

DER(1+IMAX«-IMAX) = DEVIDX(IMAX) 

CUVARI 1 ) * X 

CUVARI 2) = H, ENTHALPY - IT MUST BE H WHICH MAY BE * OR - 
CUVAR( 3 ) = Cl ( 1 ) 

CUVARI 2 + 1 MAX )=CI(IMAX) 

CUVAR(2+IMAX*1 )=EVI(1) 

CUVARI 2+IMAX+IMA X) = E VI I IMAX ) 

EXTERNAL FOFE 
C KEY INT=1 AT END OF 1ST INTERVAL 

IF I ITNEG.EQ.1.0R.IEXP.EQ.1) RETURN 
IF I M. NE.O) GO TO 10 
DO 5 1=1, IMAX 
J1 = I ♦I MAX+2 

5 CUVARI J1 )=VARIJ1 ) 

10 CONTINUE 

C INTERPOLATE FOR P ACROSS INTEGRATION INTERVAL 

DO 15 1=1, IMAX 
15 EV 1 1 I ) =CUVAR( 2*I MAX+ 1 ) 

MFTL=1 

NFTL=IZTERM— IPSI+1 

CALL FTLUP (CUVARI 1) ,PFTL, MFTL ,NFTL, VARI 1 1 PSI ) ,P I IPS I ) » 


78 


APPENDIX A - Continued 


C INTERPOLATE FOR OPOX 

CALL FTLUP (C'JVAR(1I , DPFTL ,MFT L, NF TL , VARIt I PSI ) , OPDX (IPSI II 
IF ( IBUG(IO).NE.O) GO TO 30 1 

PRINT 20. E,PFTL,P(IPSI),P<IPSI+1),CUVAR(1I,0PFTL.DP0X(IPSI> 

20 FORMAT (25H E,PFTL,P(IT ),P(IT +1/7E17.8) 

PRINT 25, ( C I ( IU), IU=1,IMAX) ,(EVI( IU >, IU=1 . IMAX) 

25 FORMAT (5E18.8) 

30 CONTINUE 

IF (KEYINT.EQ.O) GO TO 95 
SUM=0. 

C COMPUTE T VI 

DO 45 1=1, I MAX 

C USE SHOCK VALUES FOR INITIAL COMPUTATION ON EACH STREAMLINE 

SUM=SUM+CI ( 1 1 /MU 1 1 I » 

IF (EVI( D.EQ.O. I TVl(ll=0. 

IF (EVI(t).EQ.O. ) GO TO 35 

XL = I (DGENI ( I )*R*THET AKIM/I MU I( l ) *E VI ( I ) ) +1.0 > 

ALN=ALOG( XL ) 

TVI (II =THETAI( II /ALN 
35 CONTINUE 

IF ( IBUG(IOI. GT.O) GO TO 45 

PRINT 40, I.DGENII I) ,R,THETAim ,MUI (I» .EVim .XL.ALN.TVIU) 

40 FORMAT (26H I .DGENI , R, THE TA I , MU I . EV 1/ 1 5 , 6 E17. 8/ ilH XL , ALN, TVI /3 El 
18.8/1 

45 CONTINUE 

MU=1. O/SUM 
C 

U2=2.0*(HSTAG-H) 

C 

IF IU2.LT.0.I PRINT 50, H 

50 FORMAT (32H USQUARED IS NEG. IN DERSUB, H=E15.8.29H ENO STREAML 

1 INE INTEGRATION) 

IF (U2.LT.0.) IUNEG= 1 
IF (U2.LT.0. ) RETURN 
U=SQRT(U2) 

C ITERATE FOR E 

KC0DE=0 

55 CONTINUE 

K1 TR 1=0 
DELTE= 100000. 

El=. 00001 
E2=.001 
MAXI =10 
C 

CALL I TR1 (E,OELTE,FOFE»E1 ,E2, MAXI , I CODE ) 

C C5.002 LF-ITR1 

IF ( ICODE. E0.3) GO TO 75 
I.F (ICODE. EO.l) PRINT 60 

60 FORMAT (38H FOFE, MAXIMUM ITERATION EXCEEOED ***/) 

IF (ICODE. E0.2) PRINT 65 
65 FORMAT (25H FOFE, DERIVATIVES ****/> 

PRINT 70, ICODE 

70 FORMAT (13H I TR1 , I CODE= 12 ) 

IF ( ICODE. NE.l l STOP 66 
C 

C WHEN IC 0DE = 1 » TRY A NEW STARTING E. DO THIS 2 TIMES 

IF (KCODE. E0.3) STOP 66 
C 

E=EI TR1 - 

KCC»E=KC0DE+1 
GO TO 55 
75 CONT INUE 

IF ( I EXP. NE.l) GO TO 65 
E=EITR1 


STOP 66 
STOP 66 
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KC 0DE=1 +KCODE 
I E XP=0 

PRINT BO. KCODE 
80 FORMAT ( 8H KC0DE=I2) 

IF (KCODE.EO.l ) GO TO 55 
IF (KC00E.GT.2) RETURN 
E=H-0. 00001 *H 
GO TO 55 
85 CONTINUE 

RHO=PFTL/( H-E» 

T=(PFTL*MU)/IRHO*R) 

IF ( T. GT.O.OI GO TO 95 
PRINT 90, T.RHO.PFTL .H.E.MU 

90 FORMAT (14H T NEGATIVE =E 15. 8 , 2 X. 4HRH0= El 5. 8 . 2X , 3HP FTL = E15. 8 , 2X ,2 
1HH=E15.8,2X,2HE=E15.8,2X,3HMU=E15. 8) 

E=H-0. COOOO 1*H 
ITNEG=1 
RETURN 
C 

95 CONTINUE 
KE YI NT = 1 

IF ( IBUG( 9 I.NE.O I GO TO 110 
PRINT 100. IX, IPSl .RHO.T.MU.U.H 
100 FORMAT (16H I X, IP S I ,RHO , T=2 14 ,4E1 7. 8, E17. 8/ ) 

PRINT 105, E 

105 FORMAT <4H E = E18.8/( 

110 CONTINUE 

C COMPUTE DC I OX ANO OEVIOX FOR EACH SPECIE 

DO 210 I = 1 , I MA X 
C FOR EACH REACTION 

DC IOX( I |=P, 

DO 190 J = 1 , JMA X 

C SELECTED SPECIE FOR THIS REACTION 

II=MJ< J! 

IF IM.EO.O) GO TO 115 
IF ( FI (I I I.NE.O. ) GO TO 120 
115 PHI(JI=1.0 
GO TO 125 
120 CONTINUE 

TFM=EXP(THET A I ( I II/TVMIIII 
TEM1 =E XP( THETA I ( 1 1 ) / T ) 

TEM2 =F XP ( -NI ( I I ) *( THETAI ( I I I /TV I ( I I (-THETA It II l/TI I 
TEM3=EXP( THETA I( III/TVII II (-THETAI (I I) /T I 

IF (TEM3.E0.1.0I PHI ( J ) = 1. 0 
IF ( TEM3.F0.1.0I GO TO 125 

C PH I = COUPL I NG COEFFICIENT, VIBRATION ANO DISSOCIATION INTERACTION 

PHI ( J) = ( 1 1.0-TEM2) /( TEM3-1.0 )*(T EM— 1.0 )/ (TEMl-l.O) )/NI(II ) 

125 TSM4=E XPI-EJIJ l/T) 

C KJ = R AT F CONSTANT 

KJ ( J 1= A J ( J |*T**BJI J(*TEM4 
IF ( NUI J( I MAX+1, J) .EQ.O. I SJ(J) = 1.0 
IF I NUIJI IMAX+1, JI.EQ.O. ) GO TO 140 
SUM=0. 

DO 135 I SUM=1 , IMAX 
C 

SUM= SUM + AI J ( I SUM ,J )*CI ( ISUM(/MUI I I SUM) 

IF ( IBUG(7).NE.O) GO TO 135 

PRINT 130, I,J,ISUM,SUM,AIJ( ISUM.Ji.CI ( ISUMI.MUI ( ISUM) 

130 FORMAT (25H I , J , I SUM, SUM, AI J, Cl . MUI /3 1 5 , 5 E17. 8/ I 
C 

135 CONTINUE 

SJ ( J )= RHO* sum 
140 CONTINUE 

IF ( IBUG(7).NE.O) GO TO 150 
C 
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' PRINT 145. I • J.F I( II ItPHKJI ,TEM .TEM1.TEM2 .TEM3, TEM4 . AJ< J ) , 8 J( Jl «T 
1,KJ( J) .SUM.RHO.SJI J) 

145 FORMAT (26H l , J ,F I , PHI , TE M,TE Ml .TEM2/2I 5, 5E18. 8 /22H TEM3 »TEM4* AJ 
1 ,BJ,T.KJ/6E16.8/12H SUM ,RHO ,S J/3E18.8/I 
C 

150 CONTINUE 
PR0D=1.0 

DO 180 I PR0D=1 * I MAX 
IF (NUIJtIPROD.JI.EQ.O.O) GO TO 170 
TEM=RHO*CI (IPROO)/MUI( IPROD) 

IF (TEM.GE.O. ) GO T3 165 

C LET PROBLEM COMPUTE THRU DERSUB WHEN Cl NES, 

C MAKE DECISION IN MAIN ABOUT ACCEPTING NEG Cl 

ABW= AMOOI NUIJCIPROD, Jl ,2*0 » 

IF ( IBUG{8).NE.O » GO TO 160 

IF ITEM.LT. 0.) PRINT 155. I , J. I X. I PROD .TEM .RHO.MUM I PROD) ,NU I J< I PR 
10D.J ), Cl (IPROD). PROD, ABM 

155 FORMAT < 2X , 32HI , J, I X , I PROD ,TEM,RHO ,MUI ,NUI J.CI /4I3. 7E17. 8) 

160 CONTINUE 

C ABM-0.0 FOR EVEN NUIJ, ABW=1.0 FOR 000 NUIJ 

IF (ABW.EO.O.O) PR0D=PR0D* l-TEM) **NU IJ ( I PROD, J ) 

IF (ABM.EO.l.D) PROD=-1.0*PROD*(-TEM)**NUIJ< IPROD. Jl 
GO TO 170 
165 CONTINUE 

PROD=PROD*TEM**NUI J ( IPROD, J) 

C 

170 CONTINUE 

IF ( IBUGm.NS.O > GO TO 180 
IF ( IX.LT.3 ) GO TO 180 

PRINT 175, I , J , I PROD , TEM .PROD, RHO, Cl (I PROD )«NUIJ(IPROO,J) 

175 FORMAT (48H I , J , I PROD ,TEM , PROD, RHO, Cl < I PR OD ) . NU I J U PROD , J ) / 315, 6E 
117.8/) 

180 CONTINUE 

DCI JDX=PHI ( J)*KJ (J )*SJ( J )*(MUI ( I ) / < RHO*U I ) *< NUPI J( I , J )-NUI J ( I,J) I* 
IPROD 
C 

C DCIDX = RATE OF APPEARANCE AND DISAPPEARANCE OF CONCENTRATION 

C OF SPECIE 

DCIDXf I )=DC IDXI I l + DCIJDX 

C 

IF < IBUG(8).NE.O ) GO TO 190 

PRINT 185, I , J.DCI JDX.PHI ( J) ,K J< J) . S J< J) , MUI < I) . RHO. Ui NUP I J< I , J) ,N 
1UIJ( I.J), PROD, DC IDX( I ) 

185 FORMAT (25H .1 , J , DC I JDX.PHI ,K J, SJ .MUI /2 15. 5 E18. 8/ 29H RHO, U, NUPI J »N 
1UI J, PROD. DC IDX/6 E18. 8/ ) 

190 CONTINUE 

TEM=EXP(THETAI (I )/T) 

TEMl = <Fim*DGENI<I)*R*THETAI< I) )/MUI(I) 

EVIBAR ( I )=TEM1 /< TEM-1.0 ) 

IF (M.NE.O) GO TO 195 
VAR ( I*IMAX+2)»EVIBAR<I ) 

DEV I OX < 11 = 0. 

GO TO 210 
195 CONTINUE 
TAUSUM=0. 

IF ( F I ( D.EO.O.) GO TO 210 

DO 205 K=1 , IMAX 

TEM3 =E XP (-THET AI <I)/T) 

TEM4=E XP( SI GIK (K , I ) *T**( -1./3. ) I 
TEM5=(FI<I )*ALPIK<K, I ) )/PFTL 
TAUIK=TEM5*(T**BETAIK(K»II*T EM4 1 / ( 1.0-TEM3) 

TEM6=C UK)* (EVIBAR (I )-EVl( I) ) 

TAUSUM = TAUSUM*TE M6/(TAUIK*U) 

C 


APPENDIX A - Continued 


IF < IBUGI9).NE.0 ) GO TO 205 

PRINT 200, I,<,TEM3,TEM4,TEM5,TAUIK,TAUSUM,SIGIK(K.I ) .ALPIKI K, I) ,B 
1 ETAI K( K , I ) , U,E VI ( I ) , EV IBAR { I ) , TEM, TEM1 ,TEM2 , TEM6 ,TEM6 
20 C FORMAT (33H I ,K ,T EM3 , TE M4 ,TEM5 , TAUI K, TAUSUM/2 15 , 5E1 8. 8/ 33H SIGIK 

l.ALPIK ,8ETAIK,U,EVI,EVIBAR/6E17. 8/22H TEM ,TEM1 , TEM2 ,TEM6 /4E18*8 
2,022/1 

20 5 CONTINUE 
C 

C OEVIDX = EQUILIBRIUM VIBRATIONAL ENERGY 

DEVIDXU ) = T AUSUM 
OER( l^IMAX+I )=OFVIOX(I) 

210 CONTINUE 

C COMPUTE OHDX 

DHDX=DPFTL/RHO 
IF ( IBUGI 10I.NE.0) RETURN 
PRINT 215, DEV IDX< I ) 

215 FORMAT I9H DE VI OX =E 22» 14 ) 

PRINT 220, DPFTL ,DHDX 

220 FORMAT (12H DPDX, DH0X=2E18. 8, 16H END BASIC////! 

RETURN 

END 

FUNCTION fofeToum) 

C 

C CALLED BY ITR1 IN EASIC TO EVALUATE E 

C 

COMMON M* IX,IMAX, IPSI* I8UG( 15) 

COMMON P( 500 ) , DP C X ( 500),VARI( 500) 

COMMON EVIINFI25) , THETAI 1 25 ) , MU K 25 ) , F I ( 25 ) .DELHI 1 25 ) , C I INF ( 25 ) , 

1 LI (25) 

COMMON R*MUINF»DELTA»LAMSQ*OMEGSQ*OMEGA»MU*T»U 
COMMON EPSIIL(20,20),GIL(20,20) 

COMMON VAR(52),CUVAR(52),DER(51) 

C 

COMMON 2S,RS,RC. SINT, COST, PS, US.PSIS 
COMMON MJ ( 50 ) • EJ ( 50 ) , AJ ( 50 ) , BJ( 50) 

COMMON TVI(25),NI(25), DGENH25) 

COMMON NUIJ(21,30) , A I J(20,30 ) , NUPIJ(20,30) 

C 

COMMON /LAB4/ PFTL ,KITR1 , EITR1 
COMMON /LAB7/ ITNEG,IEXP 
C 

DOUBLE PRECISION VAR 
C 

DIMENSION Cl (25) 

C 

EQUIVALENCE ( CUV AR (2 ) , Cl ( 1 ) ) , (CUVAR ( 2),H) 

C 

REAL MUI * LAMSC *MU INF »MU 
C 

I F ( ITNEG.EC.l.OR.IEXP.EQ.l) RETURN 

KITR1=KITR1+1 

IF (KITR1.EQ.3) EITR1=DUK 

PORHO=H— OUM 

E=0 

EDUM=DUM 

T=PORH0*MU/R 

IF (IBUG(ll).NE.O) GO TO 10 
PRINT 5, EDUM,ML,T ,PHO,PFTL,H 

5 FORMAT ( IX , 5HEDUM=E 15.8 , 2X» 3HMU=E15. 8, 2X *2HT=E15. 8 ,2X, 4HRH0=E 15. 8, . 

1 2X,5HPFTL=E15.e,2X,2HH=E15.8) 

10 CONTINUE 
C 

DO 35 1=1, I MAX 

IF (TVI(I).EQ.O.) TEM1=(1.5*R*T)/MUI( I )+ ( F I ( I )*R*T)/MU I (I) 

IF (TVI(I).EQ.O.) GC TO 15 
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TEM=EXP( THETAI ( II/TVKI) ) 

TEM1=(1.5*R*T)/MUI( I ) + (F I ( I )*R*T)/MUI( I ) ♦( F I ( I) *R*T HET AI ( I ) ) / ( MUI ( 
II »*<TEM-1.CI » 

15 CONTINUE 
SUMG=0. 

SUMGE=0. 

C 

Lii-um 

00 25 L= 1 * L 1 1 

T EM3=— EPSI IL ( L 1 1 I/T 
IFCTEM3.LT.741.67) GO TO 42 

PRINT 41, TEH3»EPSIIL(L* 1 1» T, EDUM ,L , I ,H,MU,KITR1 

41 FORMATC 56H EXPCTEM3) , TEM3.GT.74i.67, TEM3 , EPS l IL (L , l ) , T , EOUM, L, 1= 
1 /4E17.8.2I5.13H H« MU »K I TR1=2E 1 7. 8 » I 5 ) 

1 EXP=1 

SET IEXP=1 TO TRIGGER REDUCTION IN COMPUTING INTERVAL 
AND A NEW E ESTIMATE' 

RETURN 

42 TEM2=EXP( TEM3) 

SUMG=SUMG+GIL(L, l )*TEM2 
SUMGE = SUMGE+GIL(L, I ) *EPS I IL ( L , 1 1 *TEM2 
IFCSUMGE.GT.10.E300> IBUGC1D-0 

I F t I BUG* 11) . NE . 0 ) GO TO 25 

PRINT 20, I,L,TEM2,EPSIIL(L,I),T,SUMG,GIL(L,I),SUMGE 
FORMAT (34H I , L , T EM2 ,EPSI I L , T , SUMG , G I L , SUMGE/ 2 1 5 , 6 El 8 . 8/ > 

CONTINUE 

IFCSUMGE.LT. 10. E30C) GO TO 26 
PRINT 27 

27 FORMATC 56H IN FCF E , SUMGE .GT . 10. E300 , SET IEXP=1 TO REDUCE INTERVA 
lL/> 

I BUG ( 117 = 1 
I EXP=1 
RETURN 
26 CONTINUE 

EI = TEM1+R/MUI Cl )* < SLMGE/SUMG) +DELHI C I l/MUI ( I ) 

E=E*EI*CIII> 

IF C IBUGC11 ) .NE.O > GO TO 35 

PRINT 30, I,L*IX, IPS Ij_PORHO ,E,E0UM,H,RH0,T,TEM,TEM1, SUMG, SUMGE, El, 
IE 

C FORMAT C32H I,L, IX, IPSI, PORHO, E ,EOUM,H, RH0/4I4, 6E18.8/23H T , TEM , 
1TEM1,SUMG,SUMGE/7E18.8//) 

5 CONTINUE 
FOFE=E 

IF CE.GT.H) PRINT 4C , E,H 

40 FORMATC 26H IN I TR 1-FOFE , E.GT.H, E = E15.8, 6H H=E15.8, 

1 5X, 1 7H REDUCE INTERVAL) 

IFCE.GT.H) ITNEG=1 

RETURN 
END 

FUNCTION FOFTSC CUM > 

CALLED BY ITR1 IN MAIN PROGRAM TO EVALUATE TS 

COMMON M , I > , IM AX , I PS I , I BUG 1 1 5 ) 

- COMMON PC 500) ,DPDX I 500),VARIC 500) 

COMMON EVI INFC 25 ) , TFETAI C25 ) , MU IC 25 ) ,F IC 25 ) , DELHI C 25 ) ,CI INF C 25) , 

1 LIC25) 

COMMON R.MUINF, DELTA, L AMSC , OMEGSQ, OMEGA, MU, T,U 
COMMON EPSI ILC 20 ,2C ) ,GILC20,20 ) 

COMMON VAR (52) .CUVARC52) ,0ER(51) 
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DOUBLE PRECISION VA P 
C 

DIMENSION C I (2 5 ) 

C 

EQUIVALENCE ( CUVA P ( 3 ) , C I ( 1 )) 

C 

REAL MU I » LAMS Q» ML INF f MU 
C 

E S=0 

DO 25 1 = 1 • I MAX 

IF ( M« E Q* 1 ) EVIS=E VI INF ( I ) 

IF (M.EQ.l) GO TO 5 
C 

T EM=EXP (THETAI(II/CUM) 

C 

EVIS=(R*THETAI(I) ) / ( MUI ( I) *( TEM-1 .0 )) *FI ( I) 

5 SUMG=0 

SUMGE=0 
C 

Lii=Lim 
DO 10 L= 1 • LI I 

TEM1=EXP(-EPSI I L ( L « I ) /DUM) 

SUMG=SUMG*TEMl*GIL(L,I) 

10 SUMGE=SUMGE*TEM1*GIL(L,I )*EPSIIL(L,I) 

EEIS=(R/MUI ( II )*( SURGE/ SUMG ) 

EIS=(1.5*R*DUM)/KUI( II + CFII I ) *R*DUM ) /MUI (IHEVIS + EE IS+DELHI ( I > /MUI 
1(1) 

E S=E S+EI S*C I INF( I ) 

IF ( IBUG(A).NE.O) GC TO 20 
T=DUM 

PRINT 15, IX,I,M,EVIS»TEM,T, SUMG ,SUMGE , TEM1 » EE IS , EIS, ES,CIINF(I) 

15 FORMAT ( 32H IXt I , M , EV l S ,TE M » OUM * SUMG, SUMGE/3 I 5 , 5 El 8 • 8/30H TEMl ,E 
IE IS, E IS ( I),ES,CI INF ( I )/5E 18 • 8/ 1 
C 

20 CONTINUE 

25 CONTINUE 

C 

F AC=SQRT ( OREGSO-2 **LAMSQ*( DELTA-ES ) ) 
FIN=(2.*MUINF»(DELTA-ES)<FAC)/(R*(0MEGA+FAC) ) 

FOFTS=F I N 

IF ( IBUG(4).NE.O) RETURN 

PRIN T 3 0, CM£GSQ,LA KSQ, DELTA, ES , MU INF , F AC, R , OM EGA , FIN 

30 FORMAT (44H OMEGSC , LAMSQ, DELTA, ES, MUI NF , FAC ,R , OMEGA , F IN/2( 5 El 8 . 8/ ) 
l) 

RETURN 
END 

SUBROUTINE FDPCX 

TO CCMPUTE P AND CFCX 
AND SAVE ON TAPES 9 AND 8 

COMMON M, I > , IMAX , IPS I,IBUG( 15) 

COMMON P( 5C0 ) , DP CX ( 500),VARI( 500) 

COMMON EVIINF(25),TFETAI(25),MUI(25),FI(25)»DELHI(25),CIINF(25), 

1 L I ( 2 5 ) 

COMMON R,KLINF,CELTA,LAMSC,CMEGSC,CMEGA,MU,T,U 
COMMON EPSIIL(2C,2C) , GIL ( 20 ,20 ) 

COMMON VAR(52),CUVAR(52),DER(51) 

C 

COMMON ZS,RS,RC,S I NT ,COST ,PS,US, PS IS 
CCMMCN MJ(50),EJ (50 , AJ ( 50 ) , B J( 50 ) 

COMMON TVI(25),NI(25 ),DGENI(25) 

COMMON NUIJ(21,20 ) ,A I J ( 20 , 30 ) , NUP I J ( 20 , 30 ) 

COMMON SI G IK ( 20, 20 ), ALP IK (20,20 ), BETA IK (20, 20) 

COMMON KSTAG, IMAXP2, IMAXL,N, IG, I ST AG, LPS2 , LPS 1 , IMAXP1 
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Continued 


c 

c 


5 


10 


c 

15 


C 


c 

20 


35 


40 


COMMON X,TS,ES,RHGS,HS,RHOURT,PHMAX,ClKAX,PREPSIS,PREX,PRERU 
COMMON EVIS(25)»XPST (100) 

COMMON /LAE2/ CE L X , Z STER M , I ZT ERM 

DOUBLE PRECISION VAR 

BEGIN PRESSURE DISTRIBUTION 
FOR EACH P$I FIND RANGE OF X 

00 20 IPSI=1, IZTERM 
DO 15 I X=1 , IZTERM 
1 E0F=10 

CALL RECIN ( 10 , 1 • 1 R E Cl N« X , ZS, RS , RC , COST , SI NT ) 

IF (EOF, 10) 75,5 
CONTINUE 

IF (IX. EG. IPS!) X X= X 
I E0F=1 1 

CALL RECIN ( 11 ,1 , I RECIN, I XT , T S, ES, PS ,RHOS,US, PSI S, HS ) 

CALL RECIN ( 11 ,2 , IRECIN, EV IS, 1, I MAX , 1 ) 

IF (EOF, 11) 75,10 
CONTINUE 

IF ( IXT.LT.1PSI) GC TO 15 
IF ( I XT. EG. I PS I ) PSISHK=PSIS 
IF ( lPSI.EC.IX) P(IX)=PS 
IF I I PSI • EG. IX ) GC TO 15 
P(IX) s PS+US/(RC*RS)*'(PSI SHK— PSI S ) 


CONTINUE 
REWIND 10 
REWIND 11 

CALL RECOUT (9, 1,0, XX) 

CALL RECOUT ( 9 ,2 ,0 ,P ,1 PSI , IZTERM, 1 ) 

CONTINUE 
REWIND 9 
N=1 

I E0F=9 

DO 65 IPSIM, IZTERM 
CALL RECIN ( 9 , 1 , I R EC IN , X ) 

CALL RECIN ( 9 , 2 , I R EC IN ,P , IPS I , I ZTE R M, 1 ) 

IF (EOF, 9) 75,35 
CONTINUE 
I T H I S= 0 

DO 40 I X=IPSI, IZTERM 
IF (IX.EQ.IPSI) XT=X 

IF (IX. GT.IPSI.ANC.1PSI.LT. 7) XT=XT+0ELX/5. 

IF ( IX.GT.IPSI.AND.IPSI.GT.6) XT=XT+DELX 

IF ( IX.EG.IZTERM-1) P(IX+2)=0.0 

ITHIS=ITHISn 

V AR I ( I THI S ) = XT 

CONTINUE 

I HERE=0 


C 

DO 60 IX-IPSl, IZTERM 
IF (IX. EQ. IZTERM) P(IX*l)-0. 

IF ( IX. EO. IZTERM) P(IX + 2)=0.0 
I HERE=IHERE+1 

C DIFFERENTIATE 

IF ( IX.NE.IPSI) GC TO 45 

DPDX( IX)-(-3.*P( IX )*4.*P< IX+D-PC IX+2) )/<2.*DELX) 
GO TO 50 
45 CONTINUE 

DPDX( IX)»(-P( IX-1 MPdXfll )/<2.*DELX) 

IF (IX. NE. IZTERM) GC TO 50 
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DPDX(IX)=(P(1X— 2 )-4.*P(IX-l)+3.*P( Ixn/(2.*CELX) 

50 CONTINUE 

IF (IBUG(l).NE.O) GO TO 60 
PRINT 55 * I PS I ,IX,A,DPDX( IX), P( IX) 

55 FORMAT (2I5.3E18.E) 

60 CONTINUE 

C 

CALL RECOUT ( 8 ,2 , 0 »C PDX, I PSI , IZTERM , I ) 

C 

65 CONTINUE 

. 00 70 IX=1, IZTERM 

I F ( I X.E Q. 1 > VAR I ( 1 ) -0. 

IF(IX.EQ.l) GO TO 7C 

IF (IX.LT.7) VAR I ( I X )=VARI ( IX-I )*DELX/5. 

IF (IX.GT.6) VARI(IX)=VARI( IX-1)+DELX 
7C CONTINUE 

REMIND 9 
REWIND 8 
REWIND 10 
REWIND 11 
RETURN 

75 PRINT 80, I EOF 

80 FORMAT! 7H EOF 15, 13H IN FDPDXI 

STOP 50 

C STOP 50 

END 

SUBROUTINE CHECK 

CHECK CALLED BY CALINTH TO MAKE OECISION TO ACCEPT ANSWERS 
IF ACCEPTABLE, SET ELB=0 AND RETURN 

IF NOT ACCEPTABLE MODIFY SPEC AND Cl. SET ELB=1.0 AND RETURN 

COMMON M, IX.IMAX.IPSI, IBUG(15) 

COMMON P( 5001, DPCX! 500),VARI( 500) 

COMMON EVI INF ( 25 ) ,THET AI ( 25 ) , MU I ( 25 ) , F I ( 25 ) , DELH I! 25 ) , Cl INF ( 25 ) , 

1 LI ( 25) 

COMMON R, MU INF, DELTA, LAM SC, CM EGSC, OMEGA, MU, T,U 
COMMON EPSIIL(20,2C),GIL(20,20) 

COMMON VARI52) ,CUVAP (52) ,DER( 51) 

COMMON /LAB6/ EL B , SP EC ,C J, TPR EV, HPREV , HCHECK , TCHE CK 
COMMON /LAB7/ ITNEC.IEXP 

EQUIVALENCE (CUVAR! 2),H) 

DOUBLE PRECISION VAR 

IF(IEXP.NE.l) GO TC 1 
GO TO 15 

REDUCE COMPUTING INTERVAL TO TRY TO AVOID EXP ERROR STOP 
1 CONTINUE 

IF! ITNEG.EC.O) GC TO A 

E.GT.H LEADS TC N EG T SOMETIMES LETS REOUCE INTERVAL TO TRY 

TO AVOID INSTABILITY 
GO TO 15 
A CONTINUE 

I MAXP1= I MAX+ 1 
DO 10 1=3 • I MAXP1 
C CUV AR ( 2 ) = H MAY BE NEG 

IF (CUVAR(I).LT.O.) PRINT 5, I , I X.CUVAR! I ) 

5 FORMAT ( 2AH NEG Cl, I , I X , CUVAR ( I ) = 2 15 , E 17. 8 ) 

IF (CUVAR! I) .LT.O.) GO TO 15 
10 CONTINUE 

IF (IX.LT.2) GO TC 30 

IF (IPSI. EC. 6. ANC.IX.LT. 8) GO TO 30 

IF(ABS(TPREV-T)/T .LT.TCHECK) GO TO 25 
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PRINT 11 

11 FORMAT ( 33H CHECK- AB S( TPREV-T ) /T .GE . TCHECKT ) 

15 CONTINUE 

IFISPEC.LT.1.0E-15) STOP 30 

STOP 30 

REOUCE INTERVAL 
SPEC=SPEC/4. 

C J=SPEC 
E L B= 1 • 0 

PRINT 44,SPEC,IPSI, I X, T, T PREV , H , HPRE V , I EXP, ITNEG 
44 FORMAT ( 14H RECUCEC S PEC= E15 . 8 » 10H I PSI , I X=2 1 5 , 18H T v TPREV »H» HPRE 
1V=4E16.7/12H I E XP v I TNEG=2 13) 

I TNEG=0 
I EXP=0 
RETURN 
CONTINUE 

I F ( ABS( ( HPREV-F ) /H )• GT *H CHECK J PRINT 42 
42 FORMAT ( 35H CHECK ,AE S( ( HPREV-H ) /H) . GT . HCHECKT ) 

I F ( ABS( ( HPREV-H) /H) *GT .H CHECK ) GO TO 15 
CONTINUE 

ACCEPTABLE 

TPRE V=T 
HPRE V=H 
ELB=0 • 

RETURN 

END _ __ _ 

SUBROUTINE SHOCK G ( C E LX, Z STERM , I ZTERM > 

SUBROUTINE SHOCKS C / LLEC BY MAIN 

S HOC KG SUBROLTINE MIST BE SUPPLIED BY USER. IT SHOULD DEFINE THE 
SHOCK GEOMETRY AND STORE ON TAPE 10 , X , Z S, RS ,RC ,COST , S I NT FOR EACH 
X FROM 0.0 TO X AT ZSTERM IN INCREMENTS OF CELX/5.0 TO CELX AND 
INCREMENTS OF CELX THEREAFTER 
X=DI STANCE ALONG SHOCK 
ZS=DI STANCE ALONG SHOCK AXIS OF SYMMETRY 
RS=RAD It S OF SHOCK 
RC=RADILS OF CURVATURE OF SHOCK 
COST=COS OF ANGLE CF ATTACK 
SINT=SIN CF ANGLE CF ATTACK 

USER MUST STORE THE NUMBER OF X VALUES IN I ZTERM 
USE STOP 13 FOR ERROR STOPS IN SHOCKG 


DIMENSION ZDUM(ICCO) ,RSDUM< 1000) ,RCDUM( 1000), SI DUM( 1000), 
1 CGDUMl 100C ) , XC ( 10GC ) 


NEQ A 12-16-66 
Z/L=CQSMR/L)-1.C 
L=2. 105211 CM 
NEQ2 INPUT 


EL=2. 109216 
X=0. 

DRS=. 0078125 
NRS=0 

10 NRS =NRS+1 

IF(NRS.EO.l) R SDUM ( 1 1=0. 

IF(NRS.NE.l) RSDUM ( N RS ) = RSDUM( NRS— 1 ) ♦DRS 
ZDUM(NRS) =COSH(RSCUM(NRS) 1-1.0 
XC ( NRS ) = SQRT(2.*ZtUM(NRS)+ZDUM(NRS)**2) 
CODUM (NRS ) = XC ( NRS ) /SQRT (XC( NRS ) **2+1 .0 ) 
RCDUM (NRS )= 1 .O+XC ( NRS 1**2 
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S IOUM( NRS ) = SORT ( 1 .O-CODUMC NRS )**2 ) 

RCDUM(NRS) = EL*RCCUMNRS) 

XC (NRS)=EL*XC (NRS) 

R SDUM1 NRS) = EL*RSOUMNRS) 

Z OUM( NRS )=EL*Z OUH(NRS) 

I F ( ZDUM (NRS ) .LT.ZSTERM) GO TO 10 
X=0 • 

N0X=0 
PRINT 25 

25 FORMAT ( 5X , 1HX ,16X,2HZS, 15X ,2HRS , 15X, 2HRC, 15X ,4HCQST, 13X,4HSINT/ ) 

30 NGX=NQX+ 1 

IF ( NOX. LT. 7 ) X= X + CELX/5 • 

IF (NOX.EQ.l) X=0 • 

IF (NOX.GT.6) X=X*DELX 

CALL FTLUP ( X* ZS * 1 ,NRS ,XC » ZDUM) 

CALL FTLUP < X ,R$ , 1 , NRS , XC, RSDUM) 

CALL FTLUP ( X, RC , 1 ,NRS , XC ,RC0UM) 

CALL FTLUP ( X, CO ST , 1 ,NRS , XC , COOUM ) 

CALL FTLUP < X, SINT , 1 ,NRS, XC,SI0UM) 

PRINT 35, X,ZSfRS,RC,COST,$INTtNGX 
35 FORMAT <6E17.8,I5> 

CALL RECOUT ( 10, 1 ,0, X f ZS,RS ,RC, COST, SINT ) 

IF (ZS.LT.G.O) PRINT 45 
45 FORMAT! 16H ZS IS NEGATIVE ) 

I F ( Z S • LT .0 #0 ) STOP 13 

C STOP 13 STOP 13 

IF (ZS.LT.ZSTERM) GO TO 30 
IZTERM=NOX 
ENO FILE 10 
REWIND 10 
RETURN 
END 

SUBROUTINE C AL I N T H ( V AR , DER , EL E 1 , EL E 2 , C I , SPEC , N , CU VAR , E LT , 
lCIMAXfNERRfPHMAX) 

IN THE CALINTN VERSION OF CALINT THE VARIABLE IN VAR(2) AND 

CUV AR ( 2 ) MAY BE ♦ OR -.VALUES OF OTHER DEPENDENT VARIABLES 

ARE EXPECTED TC BE PLUS. 10-69 

DOUBLE PRECISION VAP 

DIMENSION VAR(52) ,DER(51),ELE1(51) ,ELE2(51),CUVAR(52),F1(51), 
1F2(51),F3(51),CAPF1(51),CAPF2(51),CAPF3(51),P( 511 ,PH(51), 

2DELTY (51 ) , Y3 ( 51 ) • Y4 (5 1 ) , F4( 5 1 ) , Y2 ( 51 ) 

NERR=1 Cl OR N IS EQUAL TC ZERO 

NERR=2 ELE 1 LESS THAN OR EQUAL TO ELE2 


TEST INPUT 
FN=N 

TEST^C I *FN 
IF(TEST)998,997,5SE 
997 NERR= 1 
RETURN 

SS8 IF(ELEl-ELE2)999,99c f 10CC 
999 NERR=2 
RETURN 

1000 I F( SPEC ) 5 , 1 , 5 

SECTION FOR INITIALIZATION COMPUTATION OF DERIVATIVES 

1 SPEC=C I 
I CONT= 1 

2 Nl=N+l 

DO 3 1=1, N1 
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c 

3 CUVAR( I)=VAR( I) 

CALL BASIC 

RETURN WITH DERIVATIVES IN DER 

00 4 1= 1 , N 

4 F 1 ( I ) = 0ER ( I ) 


RETURN 

COMPUTE Y2,X2 

5 CUVAR<l) = VARm*CI/2. 

00 6 1 = 1 » N 

1 1=1 + 1 

Y2m-VARf I1) + CI/2.*F1CI I 
IF( I «EQ« 1 ) GO TO 6 
IF(Y2( I) >65,6,6 

6 CUVAR (ID =Y2 ( I ) 

GO TO 66 

65 SPEC=C I 
CI=CI/2. 

GO TO 5 

CALL BASIC TO EVALUATE F2 

66 CALL BASIC 

RETURN 

DO 7 1=1. N 
11 = 1+1 

F2 ( I ) = DER ( I ) 

COMPUTE Y 3 

Y3m=VAR(U) + CI/2.*F2m 
IF Cl .EQ.l) GO TC 7 
I F € Y 3 ( I) >65,7.7 

7 C UVAR ( 1 1 ) = Y 3 ( I ) 

CALL BASIC TO EV/LUATE F3 

CALL BASIC 
RETURN 

DO 10 1=1, N 
F 3 I I ) = CER ( I ) 

' COMPUTE P,PH AfvC CAP F TERMS 


*F(Y3<I)-Y2(I)) 9,8,9 

8 P(I)=0. 

GO TO 91 

9 P ( I )=-(( F3 ( I ) -F2 ( 1 ) ) /( Y3 ( I l-Y2( 1 1 1 1 
91 PH(I)=P(I)*CI 

I F (PHI I) .1 83.83,103 ... _ 

83 PH ( I ) = 0 • 

P«I)=0. 

GO TO 84 

103 IF(AB$(Y3(I)-Y2( I) I / I (ABSCY3C I H +ABSC Y2 C I ) I )/2 . > -.5E-4 ) 83. 83, 84 

84 IFCPHCI)-.l)85,85,«f 

85 CAPFi(I)=l.-PH< D/2. + CPHC I )**2)/6«-CPHCI 1**21/24. 
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C APF 2 I I > = .5-PH( IW6. + (PH( I ) **2 ) /24. - ( PH ( I ) **3 > / 120 . 
CAPF3( I )=l./6.-PHl I )/24. +(PH(I > **2 ) / 1 20 < PM I )**3)/720. 

GO TO 10 

95 CAPFlt I)=IEXP(-PH(I ))-i.)/(-PH( I)) 

CAPF2(I)=<CAPF1< I 1-1 .)/<-PH< III 
CAPF3(I )=(CAPF2( I)-.5)/(-PH( 1)1 

10 CONTINUE 

IS PH BETWEEN EL E 2 AND ELE1 

I F ( I CONT- 1)101,101,102 
102 I C0NT= IC0NT— 1 
SPEC=CI 
GO TO 17 
101 DO 11 1 = 1, N 

104 IFIAB$<PH( I) )-ELEl( I) )ll v ll»13 

11 CONTINUE 

12 SPEC=C I 
GO TO 15 

HALVE INTERVAL AND OOUBLE PH RANGE 

13 DO 96 1=1, N 
ELE1U) = ELEUI)*2. 

I FI ELEK I l-PHMAX 194,94,9 55 
94 ELE2I I ) = ELE2( I )*2 • 

GO TO 96 

955 ELEKI ) = ELE 1 ( I ) / 2 • 

96 CONTINUE 
SPEC=C I 

C I=C 1/2 • 

ICQNT=3 
GO TO 5 

RETURN TO RECOMPUTE INTERVAL 


RETURN TO RECOMPUTE INTERVAL 

15 DO 16 1 = 1, N 

IF(AB$(PH< I) )-ELE2(I ) >16,17,17 

16 CONTINUE 

DOUBLE INTERVAL 
C 1=2 • *C I 

IF(CI-CIMAX)17, 17,165 
165 C I =C I MAX 

17 CONTINUE 

COMPUTE Y 4 , X4 

DO 18 1=1, N 
11 = 1+1 

CUVARI 1 1 ) = V AR ( I1)+$PECMF3<I)*(2.*CAPF2( I ) )+Fl(I)* 
1 ( CAPF1 ( I )-2.*CAPF2(I))+F2( I ) *PH ( I ) *CAPF2 ( I ) ) 

I F ( I *EQ« 1 ) GO TO 18 
IF(CUVAR(I1))175, 18,18 
175 C I=SPEC 
C I = Cl / 2 . 

GO TO 5 

18 Y4( I ) =CUVAR ( II) 

CUVAR(1)=VAR(1)+SPEC 
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CALL BASIC TO EVALUATE F4 

CALL BASIC 
RETURN 

DO 20 1=1, N 
11 = 1*1 

F41 I ) = 0ER (II 


COMPUTE DELTA V 

DELTYU) = SPEC*(Eim*CAPFHn*(-3.*(Fim*Pm*VAR( Il» ) 
l*2.*(F2II)*P( I )*Y2(I»)+2.*(F3( I )*P( I)*Y3( I ) )-F4(I) 

2-P( I )*Y4( I ) ) *CAPF2 ( I ) *4 . * ( ( FI ( I) *P (I) * VAR 1 1 1 ) ) -( F 2 ( I ) 
3+P ( I »*Y2<m-< F3( I 1*P( I I*Y3(I» ) ♦< F4 (I) +P til * Y4( I>) ) * 
4CAPF3ID) 

COMPUTE Y*DEL T A Y 

20 CUVAR(I1) = VAR(I1)4CELTY( I) 

CALL CFECK FOR CECIS1CN TO ACCEPT OR RECOMPUTE 
INTERVAL 
CALL CHECK 
IF (ELT 121,21,23 

UPDATE Y VALUES 

21 N1=N+1 
DO 22 1=2, M 
11 = 1-1 

22 VAR(l>=VARm + DELTYl ID 
VARC 1)=VAR(1)+SPEC 

RETURN TO COMPUTE DERIVATIVES AT Y+DELTA Y 
GO TO 2 

RETURN TO RECOMPUTE INTERVAL 

23 GO TO 5 
C 

END 
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CHECK 


Test answers 


















APPENDIX A - Concluded 


The deck setup is as follows : 
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GENERAL RATE EQUATION 

In the case of a rate expression being the same for several catalytic species, it is 
desirable for the sake of convenience to incorporate all these catalysts into one reaction. 

In order to incorporate these catalysts, the set of catalytic species is treated as a general 
species designated as M. Those species denoted by M will be related to the jth reaction 
by a set of input parameters ay. Thus, if the ith species is one of the set composing M 
in the jth reaction, ay will be unity; otherwise, ay will be zero. 

If a general species M is included in the general rate expression of equation (34), 


_ K j^i /,, _ v v -rr /P c M^ I+1 ’ j 

P u \ i>j i, V iJi \ p i/ \Mm/ 


where the general species is denoted by the subscript I + 1. The mass fraction of M, 
that is, Cjyj is 


X 

= Y a y c i 


and the molecular weight is 


j. 

= ^ m iMi 


m^ being the mole fraction of the ith component of the set composing M. Hence, 


p M,i 
m; = 

P M 


p m,i = t fr pa M 


By using Dalton's law of partial pressures 


X 

= Y P M,i 
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Inclusion of equations (A2) and (A7) into equation (Al) yields a rate equation incorporating 
the general species M, 
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SPECIES, REACTIONS, AND RATES FOR CASES I TO V 
The species, reactions, and rates for Cases I to V are presented in this appendix. 

Cases I and n 

The species are N 2 , N, O 2 , O, NO, NO + , and e~ for Cases I and n. 


Number 

Reaction 

Rates 

Reference 

1 

O 2 + M — 20 + M 

1.2 x 10 21 T _3 / 2 exp^- 59 380 j 

8 

2 

20 + M -► O 2 + M 

1.0 X 10!8 T' 1 

8 

3 

NO + M — N + O + M 

5.2 x 1021 T" 3 / 2 exp(- 75 ^ 90 ) 

8 

4 

N + O + M — NO + M 

1.3 x 10 21 T" 3 / 2 

8 

5 

N 2 + M — 2N + M 

3.0 x 10 21 T" 3 / 2 exp(- 113 t 26 °] 

8 

6 

2N + M — N 2 + M 

1.67 x 10 20 T" 3 / 2 

i 

8 

7 

N + 0 2 — NO + O 

1.0 x 10 12 T 4 / 2 exp^ 3 ^°) 

8 

8 

NO + 0 — N + 0 2 

2.38 x 10 11 T^exp^- 19 13 °) 

8 

9 

N 2 + © — NO + N 

5.0 x 10 13 exp(- 38 90 °) 

8 

10 

NO + N — N 2 + O 

1.11 x 10 13 

8 

11 

N 2 + 0 2 - 2NO 

9.1 x 10 24 T' 5 / 2 exp(- - 5 900 j 

8 

12 

2NO - N 2 + 0 2 

4.79 x 10 23 T~ 3 / 2 exp^~ — 3 -^ 79 j 

8 

13 

NO+ + e" — N + O 

1.8 x 1021 T' 3 / 2 

8 

14 

N + O - NO+ + e" 

2.17 x 10 11 T°- 12 exp(- 31 858 ) 

8 
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Case III 

The species are CO2, CO, CO + , C 2 , C, O, C+, 0 + , and e" for Case m. 


mber 

Reaction 

Rates 

Reference 

1 

C0 2 + M-CO + O + M 

6.955 X 10 12 T 1 / 2 exp(- 63 286 j 

36 

2 

CO + 0 + M — C0 2 + M 

3.39 x 10 6 T 1 * 6175 

36 

3 

CO + M — C+O + M 

7.72 x 10 12 T^exp^- 128 T 925 j 

36 

4 

C + O + M — CO + M 

2.26 x 10 8 T 0,7366 exp^^-j 

36 

5 

C + O - CO + + e“ 

2.049 x 10 12 exp(- 33 30 °) 

40 

6 

CO + + e - - C + O 

6.023 x 10 16 

40 

7 

C + e‘ - C + + 2e" 

2.09 x 10 14 T 0,5 exp^- 13 ~ T 7 - 3 j 

41 

8 

C + + 2e~ ~C + e' 

9.65 x 10 21 T -0 ' 836 

41 

9 

O + e“ - 0 + + 2e“ 

2.14 x 10 14 T°- 5 exp(- 158 °°°) 

41 

10 

0 + + 2e' - O + e" 

6.05 x 10 22 T" 2 

41 

11 

C 2 + M — 2C + M 

6.74 x 10 13 T°- 5 exp(- 70 829 ) 

41 

12 

2C + M — C 2 + M 

1.61 x 10 12 T 0,743 

41 
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APPENDIX C - Concluded 
Cases IV and V 

The species are CO2, CO, C, O, A, C + , A + , and e" for Cases IV and V. 


imber 

Reaction 

Rates 

Reference 

1 

C0 2 + M-CO + O + M 

6.955 x 10 12 T^expf- 63 286 J 

36 

2 

CO + O + M — C0 2 + M 

3.39 x 10 6 T 1,6175 

36 

3 

CO + M~C+0 + M 

7.72 x 10 12 T 1 / 2 exp^- 128 925 ) 

36 

4 

C+O+M-CO + M 

2.26 x 10 8 T°- 7366 exp0^ 

36 

5 

C + e~ — C + + 2e" 

2.09 x 10 14 T°- 5 exp^- 130 713 j 

41 

6 

C + + 2e _ — C + e~ 

9.65 x 10 21 t~°- 836 

41 

7 

O + e - - 0 + + 2e" 

2.14 x 10 14 T°- 5 exp(- 158 T ° 00 j 

41 

8 

0 + + 2e' - O + e" 

6.05 x 10 22 T~ 2 

41 

9 

C + O - CO + + e- 

2.049 x 10 12 exp^- 33 300 ) 

40 

10 

CO + + e“ - C + 0 

6.023 X 10 16 

40 

11 

A + e" — A + + 2e" 

1.94 x 10M T°-5exp/- 182 83 °) 

41 

12 

A + + 2e~ — A + e” 

8.55 X 10 21 T" 1 - 075 

41 
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J R n = cosh (r /r ) - l y 
S c 3 0 S c ,0 / 


Catenary shock 


0 = 81° 

= 7.01 x 10 5 cm/sec 
Altitude = 61 km 


Ref. 8 - Body 

Ref. 8 - Streamline A 

O Present method - Body 

Present method - Streamline A 


Figure 2.- Comparison of body shape and streamline location crossing 
the shock at 0 = 81° (x/R c q = 0.16) with exact inverse method 
for catenary shock. Case I. 
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Figure 14.- Electron concentration along body surface. Case n. 
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Figure 29.- Species concentrations along 
streamlines crossing the shock at 
x = 0.4 and x=1.6. Case IV. 

















Figure 34.- Electron concentration profile Figure 35.- Ratio of forward reaction rate 
across shock layer at x = 1.6. to reverse reaction rate along inner 

Case IV. streamline (x = 0.4). Case IV. 
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Figure 38.- Temperature variation along 
streamlines crossing the shock at 
x = 0.4 and 1.6. Case V. 




Figure 39.- Species concentrations along streamlines crossing 
the shock at x = 0.4 and 1.6. Case V. 
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Figure 44.- Electron concentration profile Figure 45.- Ratio of forward reaction rate 
across the shock layer at x = 1.6. to reverse reaction rate along inner 

CaseV. streamline, (x = 0.4). Case V. 
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